# Optimization of laser stabilization via self-injection locking to a
whispering-gallery-mode microresonator: experimental study

Artem E. Shitikov  1,* Ilya I. Lykov  1 Oleg V. Benderov  2 Dmitry A.
Chermoshentsev  1,2,3 Ilya K. Gorelov  1,4 Andrey N. Danilin  1,4 Ramzil R.
Galiev  5 Nikita M. Kondratiev  5 Steevy J. Cordette  5 Alexander V. Rodin  2
Anatoly V. Masalov  1,6 Valery E. Lobanov  1 Igor A. Bilenko  1,4  1Russian
Quantum Center, 143026 Skolkovo, Russia
2Moscow Institute of Physics and Technology, 141701, Dolgoprudny, Russia
3Skolkovo Institute of Science and Technology, Moscow 143025, Russia
4Faculty of Physics, Lomonosov Moscow State University, 119991 Moscow, Russia
5Directed Energy Research Centre, Technology Innovation Institute, Abu Dhabi,
United Arab Emirates
6Lebedev Physical Institute, Russian Academy of Sciences, 119991, Moscow,
Russia
<EMAIL_ADDRESS>
###### Abstract

Self-injection locking of a diode laser to a high-quality-factor
microresonator is widely used for frequency stabilization and linewidth
narrowing. We constructed several microresonator-based laser sources with
measured instantaneous linewidths of 1 Hz and used them for investigation and
implementation of the self-injection locking effect. We studied analytically
and experimentally the dependence of the stabilization coefficient on tunable
parameters such as locking phase and coupling rate. It was shown that precise
control of the locking phase allows fine tuning of the generated frequency
from the stabilized laser diode. We also showed that it is possible for such
laser sources to realize fast continuous and linear frequency modulation by
injection current tuning inside the self-injection locking regime. We
conceptually demonstrate coherent frequency-modulated continuous wave LIDAR
over a distance of 10 km using such a microresonator-stabilized laser diode in
the frequency-chirping regime and measure velocities as low as sub-micrometer
per second in the unmodulated case. These results could be of interest for
cutting-edge technology applications such as space debris monitoring and long-
range object classification, high resolution spectroscopy and others.

††journal: oe

##  1 Introduction

The development of narrow-linewidth and highly stable lasers is one of the key
tasks of cutting-edge technologies. Such lasers provide unique opportunities
in variety of progressive areas as coherent communication [1, 2], ultrafast
optical ranging [3, 4, 5], atomic clocks [6], astronomy [7], and others.
Narrow-linewidth lasers are used in coherent Doppler LIDARs for aircraft wake
vortex measurements [8], aerosol detection [9] and remote spectroscopic
measurements [10]. Highly stable single frequency lasers are critically
important in inverse synthetic aperture LADAR systems (ISAL). ISAL is actively
used for space debris monitoring and controlling [11], long distance object
classification with spatial resolution far beyond diffraction limit of the
receiving aperture [12].

Self-injection locking (SIL) of a diode laser frequency to an eigenfrequency
of a high-Q whispering gallery mode (WGM) microresonator provides outstanding
results in laser spectral characteristics enhancement [13]. It proved to be a
simple and robust way of laser stabilization. Since it was demonstrated for
the first time [14] SIL still attracts increasing interest and still evolving.
The comprehensive theory of the SIL was developed in [15] and optimal regimes
of laser stabilization were discussed in [16]. Laser stabilization to sub-Hz
linewidth was demonstrated with crystalline microresonator in [17], and with
on-chip microresonator in [18]. Recently, simultaneous stabilization of two
diode lasers by one microresonator has been studied theoretically and
experimentally [19]. Furthermore, it was shown that SIL phenomena could be
efficiently used to demonstrate the formation of bright [20, 21, 22, 18, 23]
and dark [24, 25] microcomb solitons, ultra-low-noise photonic microwave
oscillators [26], frequency-modulated continuous wave LIDARs [27]. The
advantages of using SIL phenomena of semiconductor laser diodes with
integrated microresonators are the possibility of on-chip realization [28] of
such system that makes the technology compact and inexpensive. Lasers based on
the self-injection locking of laser diodes to integrated microresonators from
silicon nitride demonstrate outstanding performance [25, 29, 30].
Nevertheless, crystalline microresonator-based SIL lasers provides better
phase noise [31, 26] and greater opportunities for research and optimization
of the key parameters of the effect due to the flexibility unattainable for
integrated systems.

In this work we demonstrate comprehensive study of self-injection locking
phenomenon, with accurate tunability of different experimental parameters and
controllable switching between different SIL regimes. We discuss the
possibility of adaptation of the parameters of the self-injection locking
scheme for various applications. We studied the influence of the phase shift
of the backscattered wave (locking phase) and laser-to-microresonator coupling
efficiency (loading) on the performance of the self-injection locked laser
(the spectral characteristics of the resulting radiation, the stabilization
coefficient, the width of the locking range, and the resulting laser
frequency) and found out interesting opportunities for some up-to-date
applications.

For detailed experimental SIL investigation we assembled several experimental
setups with precise translation stages to vary key parameters including two
fully-packaged turn-key SIL diode lasers. The beatnote signal of two laser
diodes stabilized by high-Q MgF2 crystalline microresonators demonstrated an
instantaneous 1 Hz linewidth. Special attention was paid to the dependence of
the SIL parameters on the locking phase defined by the optical path between
laser and microresonator. Using spectrogams, we visualized experimental tuning
curves for different locking phases. Varying the locking phase, we found out
the possibility of the fast fine-tuning of laser diode generation frequency.
Also, it was revealed that, at particular value of the locking phase for
high-Q WGMs, the tuning curve splitting can be achieved.

Another parameter, that can be precisely controlled by the adjustment of the
gap between the coupler and the microresonator, is the coupling rate. We
analyzed experimental dependence of the stabilization coefficient and showed
that its maximum does not coincide with critical coupling, unlike locking
width. The experimental results are in good agreement with ours theoretical
predictions [32].

Beside frequency modulated continuous wave (FMCW) LIDARs, continuous
tunability of a laser frequency is very attractive for tunable diode laser
absorption spectroscopy (TDLAS)[33] and laser cooling. Studying tuning curves
in the self-injection locking regime, we revealed that there is an area where
the frequency changes linearly with driving current. We demonstrated that the
frequency tuning inside locking range can be realized up to dozens MHz without
dramatic linewidth degradation. We showed that this fact can be used to
realize linear frequency modulation up to 200 kHz and the amplitude of the
frequency chirp can be controlled by the locking phase. Tuning the frequency
of a laser in the self-injection locking regime in a linear way without using
an additional acousto-optic modulator seems to be an extremely attractive
opportunity for the implementation of devices such as LIDARs, spectrographs,
and optical sensors. We experimentally demonstrated the possibility of
coherent detection of a self-heterodyne signal through a delay line of 10 km
using such frequency modulation (laser chirping) in the self-injection locking
regime.

##  2 Stabilization coefficient control via loading and locking phase: theory

The stabilization coefficient $K$ is determined as the inverse derivative
(slope) of the generation frequency over the free laser frequency (laser
cavity):

$K=\left(\frac{\delta\omega_{\rm gen}}{\delta\omega_{\rm
free}}\right)^{-1}=\frac{K_{\rm free}}{K_{\rm lock}},$ (1)

where $\delta\omega_{\rm free}$ is the frequency variation in the free-running
regime, $\delta\omega_{\rm gen}$ is the generation frequency variation, which
is small in the locked regime providing high $K$. In experiment it is
convenient to present the stabilization coefficient as a ratio of the free-
running laser frequency change with the injection current $K_{\rm free}$ in
unlocked regime to the locked laser frequency change $K_{\rm lock}$. The
stabilization coefficient $K$ determines the linewidth of the locked laser
[34, 32]:

$\delta\nu_{\rm locked}=\frac{\delta\nu_{\rm free}}{K^{2}},$ (2)

where $\delta\nu_{\rm locked}$ and $\delta\nu_{\rm free}$ are the linewidths
of the locked and free-running laser, respectively.

We used the expression from [16] for the stabilization coefficient as a
function of five SIL parameters. We suppose the detuning of the laser
frequency $\zeta=2(\omega-\omega_{m})/\kappa_{m}$ from the WGM eigenfrequency
close to zero, since we consider the case of optimal stabilization [16]. Here
$\kappa_{\rm m}=\kappa_{\rm mc}+\kappa_{\rm mi}$ is the microresonator’s mode
decay rate with $\kappa_{\rm mi}$ determined by the intrinsic losses and
$\kappa_{\rm mc}$ by the coupling; $\omega$ and $\omega_{m}$ are the pump
output frequency and eigenfrequency of the microresonator. In this case, one
can obtain the following expression for the stabilization coefficient:

$K=1+4\frac{\bar{\kappa}_{\rm do}}{\kappa_{\rm
mi}}\eta(1-\eta)\beta\frac{(2+\frac{\kappa_{\rm
m}\tau_{s}}{2}(\beta^{2}+1))\text{cos}\psi}{(\beta^{2}+1)^{2}},$ (3)

where $\bar{\kappa}_{\rm do}$ is the laser effective output beam coupling
rate. The corresponding values of the loaded and intrinsic quality factors are
defined as $Q_{\rm m}=\omega/\kappa_{\rm m}$, $Q_{\rm int}=\omega/\kappa_{\rm
mi}$ and $\eta=\kappa_{\rm mc}/\kappa_{\rm m}$ is the coupling coefficient
characterising microresonator loading, $\beta=2\gamma/\kappa_{\rm m}$ is the
normalized mode-splitting coefficient [35], $\gamma$ is the forward-backward
wave coupling rate, $\tau_{s}$ is the feedback round-trip time.

Figure 1: Sketches of the experimental setups: (a) General optical scheme,
which includes block "SIL laser" consisting of a laser diode without isolator
exciting WGMs through a coupling prism, and RF source that varies laser
injection current. That "SIL laser" block presents in the following setups.
(b) Heterodyning of SIL laser with reference laser to figure out the
dependence from the locking phase. (c) Self-heterodyning linewidth measurement
scheme with delay line. (d) Heterodyning of two different SIL lasers. (e)
Reflector velocity measurement setup with SIL laser. Here ESA is the
electrical spectrum analyzer, OSA is the optical spectrum analyzer, OSC is the
oscilloscope, PD is the $\rm 5\ GHz$ bandwidth photodetector, RF source is the
source for LD injection current modulation.

Let us take the following equation from [36] for the critical coupling regime
($\kappa_{\rm mi}=\kappa_{\rm mc}$):

$\bar{\kappa}_{\rm do}=\frac{4}{3\sqrt{3}}\delta\omega_{\rm
crit}\frac{\kappa_{\rm mi}}{\gamma},$ (4)

where $\delta\omega_{\rm crit}$ is the width of the locking range of the
stabilized laser in the critical coupling regime. This equation can be
substituted in (3). After that considering the limit of small $\beta\leq 0.1$,
that is exact true for high-Q crystalline microresonators in the most cases,
one can simplify (3):

$K=1+\frac{64}{3\sqrt{3}}\delta\omega_{\rm crit}\frac{\kappa_{\rm
mc}\kappa_{\rm mi}}{(\kappa_{\rm mi}+\kappa_{\rm mc})^{3}}\text{cos}\psi.$ (5)

If critical coupling is assumed, then (5) will take the following form:

$K=1+\frac{8}{3\sqrt{3}}\frac{\delta\omega_{\rm crit}}{\kappa_{\rm
mi}}\text{cos}\psi.$ (6)

That equation directly bounds stabilization coefficient and locking phase, the
parameters that can be measured in experiment.

In Eq. (5) $\kappa_{\rm mc}$ depends on the gap between the coupling element
(prism) and the microresonator as follows [37]:

$\kappa_{\rm mc}\approx\frac{\omega}{2}\dfrac{{\rm
exp}{\\{-2kd\sqrt{n^{2}-1}\\}}}{(\dfrac{n^{2}-1}{n}ak)^{3/2}}\times\begin{cases}\frac{1}{\sqrt{\pi/(1+\sqrt{n^{2}-1})}}\text{,}\;p=0\\\
1/(\sqrt{2p}\pi)\;\;\;\;\;\;\;\text{,}\;p>0\end{cases}$ (7)

where $a$ is the radius of a microresonator, $\omega$ is the optical pump
frequency, $k=\omega/c$ is the wavenumber, $p\in\mathbb{N}$ is the WGM
vertical index, $d$ is the prism-microresonator gap; this equation is valid in
case of close values of refractive index (n) of microresonator and coupling
element. It is worth noting that this equation was obtained for the case of
the spherical microresonator, but it was shown that for the spheroid
microresonator the given equation can be used with great accuracy [38, 36].
Substituting Eq. (7) into Eq. (5), we then obtain the stabilization
coefficient as a function of $d$.

Thus, one can see from obtained analytical formulas that spectral
characteristics of the laser sources stabilized via SIL effect depends on the
locking phase value and on loading defined by the prism-microresonator gap.
Further, these dependencies will be checked experimentally.

##  3 Experimental setups

A 1550 nm distributed feedback (DFB) laser with power up to $1.5$ mW is used
to excite the WGMs in bulk crystalline MgF2 microresonator. The injection
current of the laser diode can be controlled by an external generator. A
schematic sketch of the experimental setup is shown in Fig. 1(a). The laser
with a lens is installed on a high-precision translation stage with PZT, which
allows us to control the distance between the laser and the resonator cavity
with an accuracy of tens of nanometers. The resonator is also set to a
translation stage with PZT with a step of 27 nm to control the distance
between the resonator and the coupling prism. Laser injection current is
controlled with RF generator. These elements (laser diode, RF generator,
coupling prism, microresonator and lenses) compose "SIL laser" block. This
block is also used as a component of the following setups. After the coupling
element, the laser beam is coupled into a single-mode optical fiber and
proceeded to measurement equipment.

First, we characterized the self-injection locked laser diode by beating its
output signal with the reference laser (Koheras Adjustik with instantaneous
linewidth of less than 100 Hz) providing heterodyne scheme [see Fig. 1(b)].
The wavelength of the reference laser is matched to the wavelength of the
laser under test to achieve microwave beatnote. The result beatnote signal was
analyzed with a 5 GHz broadband detector, which signal was observed at an
electrical spectrum analyzer (ESA) and an oscilloscope. This setup was used to
visualize tuning curves for different locking phases.

We make noise characterization of the self-injection locked laser diodes. In
case of the absence of laser much narrower than the laser under study, the
optimal way is to heterodyne it with an equivalent laser. Fitting the obtained
beatnote with a Voigt profile one may estimate the laser linewidth. We
measured phase noise spectral density of the beatnote of two SIL laser diodes
by heterodying it [see Fig. 1(d)]. Signal was measured with electrical
spectrum analyzer with an option of quadrature analysis through detector with
bandwidth of 45 GHz.

Fully characterized lasers are supposed to be ideal candidates for some proof-
of-concept experiments. One of it is to measure a delay line length using
linearly modulated laser in SIL regime [see Fig. 1(c)]. The laser radiation is
divided into two parts, the first goes through an optical fiber about 40 cm
long with a fiber polarization controller in it, the second along a 10 km long
delay line, then the 2 radiation parts are then summed through another fiber
beam splitter and fed into a detector. The signal from the fiber detector can
be examined with an ESA and an oscilloscope. We were able to use a frequency
modulated continuous wave scheme by modulating the laser injection current and
consequently modulate output frequency inside SIL regime. The frequency
modulation parameters can be manipulated by varying locking phase.

Another application is a velocity sensor based on a Doppler frequency shift
[see Fig. 1(e)]. SIL laser diode beam is amplified with an optical amplifier
up to 30 mW. Then light through a circulator goes to a free space where
reflects from the target with diffusion surface. The target is mounted on the
piezoelectric stage which is continuously moving with low constant velocity of
several $\mu$m per second. That velocity can be measured by detecting a
Doppler-shifted frequency.

##  4 Locking phase tuning

We performed experimental investigation of the modification of SIL dynamics
and stabilized laser parameters upon locking phase tuning.

One may change locking phase varying the distance between the laser diode and
microresonator. The injection current of the diode laser was swept to obtain
generation frequency in a vicinity of the WGM eigenmode, so that self-
injection locking took place. Calculating spectrograms of the beatnote signal
allows to visualize a tuning curve for different locking phases and provides
detailed information about the self-injection locking dynamics. A spectrogram
is a two-dimensional representation of the dependence of the power spectral
density of the signal on time, with the signal intensity represented in color.
The spectrograms are calculated using the window Fourier transform with the
Blackman window.

Figure 2: SIL dynamics at locking phase $\psi=0$. The transmission signal is
presented in the top graph (a). High-order mode appeared at $\approx$ 11 ms.
The spectrogram of the beatnote signal with the reference laser is presented
in the middle (b). Enlarged areas of the spectrograms corresponding to locking
ranges are presented in the panels (c, d). The modeling results are presented
in the bottom panels (e, f). Analytical stationary solutions are plotted with
blue dashed line and dynamic modeling with solid green line.

The spectrograms of the transmission spectrum for optimal locking phase
($\psi=0$) and non-optimal locking phase ($\psi=\pi$) are presented in Figs. 2
and 3 as the most demonstrative ones to analyse the SIL changes during locking
phase changes. Laser diode current is monotonously increased from 0 to 8 ms
and decreased from 8 to 16 ms, which corresponds to frequency decreasing and
frequency increasing consequently. If there were no SIL the laser frequency
would pass through zero frequency, where it is equal to heterodyne frequency,
and then increase until the turning point at 8 ms. That is because we
calculate the absolute value of the frequency with Fourier transform. In Fig.
3 the tuning range is larger than in Fig. 2. The SIL regions are marked with
light green stripes above the panels and circumscribe with dashed lines.
Transmission spectra are presented in the top panels (a) of Figs. 2 and 3. The
consequent spectrograms are presented in the middle plots (b) and enlarged
areas of the SIL regime are presented in the bottom panels (c), (d). Almost
equal widths of the locked regions for the laser frequency increasing and
decreasing is supposed to be an indicator of $\psi=0$ [see Fig. 2]. A high-
order mode is also excited at 11 ms in the (a) panel. It is pronounced in one
direction of the frequency scan because of the SIL which takes place
asymmetrically considering microresonator eigenfrequency, see Fig. 2 (e, f).
It can be noticed in the spectrogram in the panel (b) at 11 ms as a small
plateau. The case of $\psi=\pi$ is also symmetric regarding the current sweep
however, the frequency change in the locking range is a lot higher, indicating
lower stabilization coefficient[32]. Furthermore, the splitting of the tuning
curve can occur in this regime [see Fig. 3(b,c,d)]. The position of the jump
is utterly sensitive to the locking phase. The value of the jump depends on
the quality factor, and it disappears in an overcoupled regime when
$\kappa_{\rm mc}$ becomes dominant, but persists in the undercoupled regime
when the quality factor is still high. This phenomenon can be explained as
follows: due to the high Q factor of the microresonator, even weak amplitudes
of the incoming waves can be resonantly amplified and backscattered (with a
rate about 1/150) to stabilize the laser emission. In other words the feedback
from the edge of laser emission line is resonantly reflected and initiates the
transition. So the finite laser emission linewidth causes the spontaneous
transition of the laser emission frequency to the branch with higher
stabilization coefficient [see Figs. 2,3(e, f)] before the turning point. We
performed a dynamical modeling of the self-injection locking[32] and found
that such spontaneous transition is possible after certain microresonator
quality factor value [see green lines Figs. 2,3(e, f)]. Similar spontaneous
locking was also seen in nonlinear regime in [22]. We note that laser rate
equations automatically generate nonzero linewidth in the modelling due to
numerical noise. As can be seen in Fig. 3, the experimental observations (c,d)
are in full agreement with the theoretical simulations of SIL dynamics (e, f)
and the corresponding splittings for forward and backward scan in the case
$\psi=\pi$ are present.

Figure 3: SIL dynamics at locking phase $\psi=\pi$. The transmission signal is
presented in the top graph (a). The spectrogram of the beatnote signal with
the reference laser is presented in the middle (b). Enlarged areas of the
spectrograms corresponding to locking ranges are presented in the panels (c,
d). The modeling results are presented in the bottom panels (e, f). In the
inset enlarged areas of the transmission curves are presented. Analytical
stationary solutions are plotted with blue dashed line and dynamic modeling
with solid green line.

One may use spectrograms to measure the stabilization coefficient according to
(1) as a ratio of the free-running laser frequency change with the injection
current $K_{\rm free}$ [e.g. line slope in unlocked region in Fig. 2(b)] and a
locked laser frequency change $K_{\rm lock}$ [e.g. line slope in Fig.
2(c),(d)]. The values for $K_{\rm free}$ and $K_{\rm lock}$ for different
tunable parameters were obtained from spectrograms. We make linear
approximation of the frequency tuning of the laser in unlocked regime for more
than hundred independent measurements and calculate average value of the
$K_{\rm free}\,=\,(44\pm 4)$ MHz/ms. Then, we approximated frequency
dependence in locked states with fifth order polynomials and find minimum of
the frequency slope that corresponds to maximum $K_{\rm lock}$, see locked
areas in Fig. 2, 3(c, d). The ratio of the obtained values of $K_{\rm free}$
and $K_{\rm lock}$ allows us to determine the stabilization coefficient $K$,
as defined in equation (1). Thus, an experimental dependence of the
stabilization coefficient on the gap (coupling) and locking phase are
obtained. The linewidth of a free-running laser is approximately 2 MHz. The
squared stabilization coefficient maximum is $(57\pm 5)\cdot 10^{3}$, so
according to Eq. (2) we expect locked linewidth of $\approx$ 40 Hz, additional
information about calculating stabilization coefficient from Eq. (3) is
presented in supplementary materials.

To analyze how the stabilization coefficient of the SIL changes with the
locking phase variation [see Fig. 4] we measure beatnotes for the spectrograms
for different locking phases in the laser diode displacement range of $0.8\
\mu$m. It corresponds to phase rotation for more than $2\pi$ since it changes
the optical path for more than $\lambda$. The stabilization coefficient
dependence on the locking phase is calculated from the spectrograms according
to Eq. (1) and compared with values obtained using (6) with $\delta\omega_{\rm
crit}$ taken from experimental data. The experimental results are in a good
agreement with theoretical curve based on Eq. (6) shown in Fig. 4(a).

Figure 4: (a) The dependence of the stabilization coefficient on the locking
phase. The theoretical curve is calculated using Eq. (6) and data presented as
dots are obtained from the spectrograms using Eq. (1) ($\delta\omega_{\rm
crit}/(2\pi)=610\pm 20$ MHz and $\kappa_{\rm mi}/(2\pi)=3.73\pm 0.22$ MHz,
which corresponds to the microresonator intrinsic quality factor $Q_{\rm
int}=(5.1\pm 0.3)\cdot 10^{7}$). (b) The dependence of the locking range width
on the locking phase.

The locking range is also dependent on the locking phase [see Fig. 4(b)]. We
measure locking range that is sum of the locking ranges during frequency
increasing $\delta\omega_{\rm right}$ and decreasing $\delta\omega_{\rm
left}$: $\delta\omega_{\rm lock}=\delta\omega_{\rm right}+\delta\omega_{\rm
left}$. It can be seen that the period of the stabilization coefficient $K$
and locking width from the locking phase coincide. The maximum of K is
localized near the $\psi=0$ phase, where the widths of the locking ranges at
both directions of the frequency scan become equal. The nonsmoothness of the
curves can come from the fact that the locking width is not determine well due
to spontaneous locking.

##  5 Coupling optimization

Another important controllable parameter of the SIL is the coupling rate. The
simple relationships that describe the dependence of the SIL laser spectral
characteristics on the distance between the coupling element and the WGM
microresonator is obtained in Section II (5), (7). We also developed an
original setup [see Fig. 1(a)] and performed experimental verification of the
obtained theoretical results.

We use spectrograms to measure dependencies of the stabilization coefficient
and locking width on a gap between the coupling prism and a microresonator.
Locked areas expand until critical coupling and get narrower during
overcoupling. We measure locking width and stabilization coefficient for the
different values of the gap with 27 nm step. The dependence of SIL parameters
on the gap are measured for the optimal locking phase near $\psi=0$, and its
value remains constant during loading as the gap is not included into optical
path.

Figure 5: The dependence of the locking range (a) and the stabilization
coefficient (b) on the gap between the coupling element and the
microresonator. The experimental data are plotted with dots and theoretical
estimations according method from [36] and Eq. (5) are plotted with solid
lines. Approximation parameters are close and overlapped with error deviations
$\kappa_{\rm mi}/(2\pi)=(3.73\pm 0.22)\cdot 10^{6}$ Hz for the locking width
and $\kappa_{\rm mi}/(2\pi)=(3.82\pm 0.09)\cdot 10^{6}$ Hz for the
stabilization coefficient. One may see the extremum of the stabilization
coefficient shifted to the undercoupled regime, which is exact consequence of
Eq. (5).

The dependencies of the locking range width and stabilization coefficient on
the gap between coupler and microresonator are presented in Fig. 5 in the top
and bottom panels consequently. We approximate the locking range dependence on
the gap using the method described in [36]. The approximation parameters were
determined as follows: vertical index of the mode $p=2$, $\kappa_{\rm
mi}/(2\pi)=3.73\pm 0.22$ MHz $(Q_{\rm int}=(0.51\pm 0.06)\cdot 10^{8})$ and
locking range at critical coupling $\delta\omega_{\rm crit}/(2\pi)=610\pm 20$
MHz. The maximum of the locking range corresponds to the critical coupling.

The experimental values of the stabilization coefficient obtained from the
spectrograms were approximated by a theoretical curve (Eq. (5)). The
approximation parameters were determined as follows: vertical index of the
mode $p=2$, $\kappa_{\rm mi}/(2\pi)=3.82\pm 0.09$ MHz $(Q_{\rm int}=(0.50\pm
0.01)\cdot 10^{8})$. The value of $\kappa_{\rm mi}$ was used for the
calculations of the separately measured stabilization coefficient dependence
on the locking phase in Fig. 4(a) using Eq. (6). The locking range presented
in Fig. 4(b) at the optimal phase coincides with the locking range at critical
coupling in Fig. 5(a) that indicates optimal phase during the experiment with
coupling variation. Finally, the maximum stabilization coefficient from Fig.
4(a) 239 $\pm$ 15, and maximum $K$ from Fig. 5(b) exceeds 270, but at the
critical coupling it is from 215 to 256, that is perfect agreement. Such a
close relationship between the parameters obtained independently in various
experiments indicates a high reliability of the theory.

It is worth discussing that maximum of the stabilization coefficient does not
coincide with critical coupling, that is completely according to theoretical
estimations. It was shown in [16] that the optimal coupling condition for low
$\beta$ tends to $\eta=1/3$. The value of the stabilization coefficient is
decreased by 20% at critical coupling and since the laser linewidth is inverse
proportional to the stabilization coefficient squared it seems even more
important. Also, consequently, any technical deviations, for example, of the
locking phase or injection current, will be suppressed more effectively in the
undercoupled regime than in critically coupled one that will influence not
only the instantaneous linewidth but the relatively long-term stability. In
the case of the perfect mode matching the output from the prism for critical
coupling is equal to zero. So, the obtained optimal value is appropriate for
applications where the uncoupled light is used.

##  6 SIL laser characterization

Figure 6: Single sideband spectral density of the relative phase noise of the
two equivalent SIL laser diodes is presented in the top panel. The beatnote
spectrum is presented below. The enlarged area of the peak with Voigt fit is
presented in the inset, Lorentzian part is 25 Hz and Gaussian part is 540 Hz.
Linewidth measurements performed with 100-Hz resolution bandwidth.

An experiment was carried out to measure the beatnote signal of two SIL lasers
stabilized with high-Q WGM microresonators. The coupling was set up close to
the critical. The injection current was set up close to the laser generation
threshold in both cases to avoid nonlinearities and reduce thermal frequency
deviations. Varying injection current and laser diode temperature we found
WGMs which eigenfrequencies were different from each other by 3 GHz and both
modes simultaneously provided effective stabilization. The loaded quality
factors of the modes were approximately $10^{8}$. Single sideband spectral
density (SSB) of the beatnote phase noise is presented in the top panel of
Fig. 6. At the frequency offset of $5\cdot 10^{4}$ Hz phase noise level
corresponds to instantaneous linewidth about 1 Hz. At higher frequency offsets
the noise level of the measuring equipment become dominant, see the plateau
from $10^{5}$ Hz.

The beatnote signal is presented in the bottom panel of Fig. 6. The
approximation of the beatnote is presented in the inset. Linewidth
measurements performed with 100-Hz resolution bandwidth. To analyze the beat
data, we use an approximation with the Voigt profile (black curves in the
panels). The Voigt profile allows us to estimate of the contribution of the
white noise (responsible for the Lorentzian linewidth) and the flicker noise
impact (Gaussian linewidth). The Gaussian part of the approximation is equal
to 540 Hz and Lorentzian part is 25 Hz.

##  7 Linear frequency modulation

Figure 7: Spectrum of the signal of the self-beatnote with the delayed light
during coherent optical ranging. The laser was modulated inside SIL regime
with the period of 1 ms. The peak, locateded at 150 kHz, corresponds to
amplitude of the frequency modulation of 1.5 MHz given the 10 km delay line,
that is in a good agreement with frequency changing inside SIL [see Fig.
2(c,d)].

For many important purposes it is necessary to realize linear frequency
tuning. The ability to manipulate the frequency linearly by controlling the
injection current of a laser diode was shown for conventional diode lasers.
However, due to instabilities, special tricks or additional equipment are
required to linearize the tuning, for example, machine learning can be used
for this purpose to compensate the nonlinear frequency response for a
particular laser diode by restructuring the injection current [39].

In Fig. 2(c,d) one can find the areas where the laser frequency changes
linearly with the injection current in the SIL regime. Moreover, these
sections are simultaneously the regions of the highest stabilization
coefficient. The deviation of the frequency changes in this case is 1-10 MHz
and depends on the position on the tuning curve, locking phase, and coupling.

We observed linear frequency modulation with frequencies from 150 Hz to 200
kHz. The resonator was critically coupled, and the locking phase was close to
optimal one. The maximum amplitude of the injection current modulation depends
on the locking range width, the wide range is the higher injection current
modulation inside SIL can be achieved. The injection current modulation is
applied by means of RF-generator.

We conducted experiment to demonstrate conceptually a FMCW LIDAR in the SIL
regime. The laser frequency was modulated linearly inside the locking range
[e.g. with amplitude 1.5 MHz, see Fig. 2(c,d)] with a period of 1 ms. Note
that the frequency increases and decreases with the rate of 3 MHz/ms for that
amplitude during the period. The laser beam was split into two parts, one was
propagating through a 10 km delay line and the other – only through a
polarization controller. The spectrum of the self-heterodyne signal is shown
in Fig. 7. To calculate the spectrum, we took the signal sections minus 100
microseconds near the turning point which appear when the frequency tuning
changes direction. The length of the optical fiber is $l=10$ km which
corresponds to the time delay of $\tau=ln_{\rm eff}/c\approx 50\ \mu$s, where
$n_{\rm eff}$ is effective refractive index of the optical fiber. The
frequency difference between the beams caused by the time delay is expected to
be 150 kHz as a delay multiplied on the frequency changing rate, which is in
good correspondence with the measured spectrum. If we will turn off the self-
injection locking the self-heterodyne signal will disappear.

##  8 Velocity measurements

Diode lasers are widely available efficient sources of laser radiation,
however, for measuring velocities of the order of $\mu$m/s, the use of these
lasers is limited due to the linewidth of the order of 1 MHz, as discussed,
for example, in [40], as well as in [41]. SIL lasers seem to be ideal
candidates for measuring low velocities, e.g. via Doppler shift measurement,
due to achievable ultra-narrow linewidths.

The experimental setup for the measurement of the Doppler frequency shift of
the reflected signal of the SIL laser is presented in Fig. 1(e).

Figure 8: Reflected signal from the moving target. (a) The signal in case of
velocity of the target of 1.4 $\mu$m/s over 20 s measurement time. (b) Fourier
transform of signal from (a). The maximum is observed at 1.67 Hz which is
close to expected 1.75 Hz. (c) The sector of the signal fitted with 1.75 Hz.
(d) The signal in case of velocity of the target 0.27 $\mu$m/s over 100 s. The
signal is noisier than for 1.4 $\mu$m/s. (e) Fourier transform of signal from
(d). The maximum is observed at 0.29 Hz which is close to expected 0.35 Hz.

Diffuse reflective plate was used as a target. It was placed in the focus of a
converging lens. The plate was fixed on the translation stage with a PZT. Ramp
voltage with amplitude of 4 V that corresponds to 6.8 $\mu$m displacement and
frequency of tens of mHz was applied to the PZT. As can be seen from the
signal, in the vicinity of the turning point, the frequency of the sinusoidal
changes [see Fig. 8(a)]. Fourier transform of this signal gives us maximum at
1.67 Hz [see Fig. 8(b)]. The frequency of the signal is close to expected
(1.75 Hz for velocity of 1.4 $\mu$m/s, if we extract the part corresponding to
the uniform motion of the surface between the turning points. We extracted the
sector from (a) and fitted it with sinusoidal function with 1.75 Hz [see Fig.
8(c)]. It is obvious that if we turn off SIL by receding the microresonator
the frequency cannot be measured in that way due to high level of phase noise.

The lowest velocity we observed was 0.27 $\mu$m/s [see the signal in Fig.
8(d)]. The maximum of the Fourier transform was located at 0.29 Hz instead of
0.35 Hz expected, that corresponds to 0.22 $\mu$m/s instead of 0.27 $\mu$m/s.
This variance can be caused by roughness of the target which can be pronounced
at this level of measured velocities. Nevertheless, demonstrated accuracy
conceptually proves the effectiveness of SIL lasers for such sort of
measurements.

##  9 Conclusion

We investigate in detail such aspects of the self-injection locking phenomena
as dependence of the spectral characteristics on the locking phase and
loading. We have verified experimentally the dependence of the stabilization
coefficient on the locking phase, found the optimum and showed that the
resulting linewidth remains narrow in a wide range of phases. We have shown
that one can vary the locking phase within the limits defined by the desired
linewidth and thus tune the laser frequency in a range of several MHz. Due to
the purely optical essence of SIL the tuning can be very fast and, in our
case, limited by the PZT drive. The phenomenon of spontaneous locking was
observed and its bad influence on the locking at $\pi$-phase was revealed. We
demonstrate that the locking phase does not change with loading. It was
demonstrated experimentally that the maximum of the stabilization coefficient
appears in the undercoupled regime and does not coincide with critical
coupling. The SIL lasers used in experiments demonstrate superior stability,
its phase noise single sideband spectral density was measured by heterodyning
of two equivalent laser sources and 1 Hz instantaneous linewidth was
demonstrated analyzing phase noise spectral density and 25 Hz linewidth by
beatnote approximation with Lorentzian profile. As an example of possible
applications, we demonstrated linear frequency modulation by the injection
current manipulation inside SIL regime. The repetition rate was up to 200 kHz.
This modulation was used to measure the length of a $\simeq 10$ km delay line.
The SIL laser source was also used for sub-$\mu$m velocity measurements. These
experiments demonstrate new ways of the SIL scheme optimization and
implementation in a fashion of the perspective of application SIL lasers for
long distance measurements, space debris monitoring and controlling, distant
object classification and other applications.

Funding RQC team was supported by the Russian Foundation for Basic Research
(Project 19-29-06104mk).

Acknowledgments AES and VEL acknowledge the Foundation for the Advancement of
Theoretical Physics and Mathematics “BASIS” for personal support.

Disclosures The authors declare no conflicts of interest.

Data availability Data underlying the results presented in this paper are not
publicly available at this time but may be obtained from the authors upon
reasonable request.

Supplemental document See Supplement 1 for supporting content.

## References

  * [1] A. Fülöp, M. Mazur, A. Lorences-Riesgo, Ó. B. Helgason, P.-H. Wang, Y. Xuan, D. E. Leaird, M. Qi, P. A. Andrekson, A. M. Weiner, and V. Torres-Company, “High-order coherent communications using mode-locked dark-pulse Kerr combs from microresonators,” Nature Communications 9, 1598 (2018).
  * [2] P. Marin-Palomo, J. N. Kemal, M. Karpov, A. Kordts, J. Pfeifle, M. H. P. Pfeiffer, P. Trocha, S. Wolf, V. Brasch, M. H. Anderson, R. Rosenberger, K. Vijayan, W. Freude, T. J. Kippenberg, and C. Koos, “Microresonator-based solitons for massively parallel coherent optical communications,” Nature 546, 274–279 (2017).
  * [3] P. Trocha, M. Karpov, D. Ganin, M. Pfeiffer, A. Kordts, S. Wolf, J. Krockenberger, P. Marin-Palomo, C. Weimann, S. Randel, W. Freude, T. Kippenberg, and C. Koos, “Ultrafast optical ranging using microresonator soliton frequency combs,” Science 359, 887–891 (2018).
  * [4] M.-G. Suh and K. Vahala, “Soliton microcomb range measurement,” Science 359, 884–887 (2018).
  * [5] J. Riemensberger, A. Lukashchuk, M. Karpov, W. Weng, E. Lucas, J. Liu, and T. Kippenberg, “Massively parallel coherent laser ranging using a soliton microcomb,” Nature 581, 164–170 (2020).
  * [6] Z. L. Newman, V. Maurice, T. Drake, J. R. Stone, T. C. Briles, D. T. Spencer, C. Fredrick, Q. Li, D. Westly, B. R. Ilic, B. Shen, M.-G. Suh, K. Y. Yang, C. Johnson, D. M. S. Johnson, L. Hollberg, K. J. Vahala, K. Srinivasan, S. A. Diddams, J. Kitching, S. B. Papp, and M. T. Hummon, “Architecture for the photonic integration of an optical atomic clock,” Optica 6, 680 (2019).
  * [7] M.-G. Suh, X. Yi, Y.-H. Lai, S. Leifer, I. Grudinin, G. Vasisht, E. Martin, M. Fitzgerald, G. Doppmann, J. Wang, D. Mawet, S. Papp, S. Diddams, C. Beichman, and K. Vahala, “Searching for exoplanets using a microresonator astrocomb,” Nature Photonics 13, 25–30 (2019).
  * [8] S. Wu, X. Zhai, and B. Liu, “Aircraft wake vortex and turbulence measurement under near-ground effect using coherent Doppler lidar,” Optics Express 27, 1142–1163 (2019).
  * [9] G. Dai, X. Wang, K. Sun, S. Wu, X. Song, R. Li, J. Yin, and X. Wang, “Calibration and retrieval of aerosol optical properties measured with coherent Doppler lidar,” Journal of Atmospheric and Oceanic Technology 38, 1035–1045 (2021).
  * [10] K. Wang, C. Gao, Z. Lin, Q. Wang, M. Gao, S. Huang, and C. Chen, “1645 nm coherent Doppler wind lidar with a single-frequency Er: YAG laser,” Optics Express 28, 14694–14704 (2020).
  * [11] D. Mo, N. Wang, G. Li, R. Wang, and Y. Wu, “3-D inverse synthetic aperture ladar imaging and scaling of space debris based on the fractional Fourier transform,” IEEE Geoscience and Remote Sensing Letters 16, 236–240 (2018).
  * [12] N. Wang, R. Wang, D. Mo, G. Li, K. Zhang, and Y. Wu, “Inverse synthetic aperture LADAR demonstration: system structure, imaging processing, and experiment result,” Applied optics 57, 230–236 (2018).
  * [13] R. Galiev, N. Pavlov, N. Kondratiev, S. Koptyaev, V. Lobanov, A. Voloshin, A. Gorodnitskiy, and M. Gorodetsky, “Spectrum collapse, narrow linewidth, and Bogatov effect in diode lasers locked to high-Q optical microresonators,” Optics express 26, 30509–30522 (2018).
  * [14] V. Vassiliev, V. Velichansky, V. Ilchenko, M. Gorodetsky, L. Hollberg, and A. Yarovitsky, “Narrow-line-width diode laser with a high-Q microsphere resonator,” Optics Communications 158, 305 – 312 (1998).
  * [15] N. M. Kondratiev, V. E. Lobanov, A. V. Cherenkov, A. S. Voloshin, N. G. Pavlov, S. Koptyaev, and M. L. Gorodetsky, “Self-injection locking of a laser diode to a high-Q WGM microresonator,” Optics Express 25, 28167 (2017).
  * [16] R. Galiev, N. Kondratiev, V. Lobanov, A. Matsko, and I. Bilenko, “Optimization of laser stabilization via self-injection locking to a whispering-gallery-mode microresonator,” Phys. Rev. Applied 14, 014036 (2020).
  * [17] W. Liang, V. Ilchenko, D. Eliyahu, A. Savchenkov, A. Matsko, D. Seidel, and L. Maleki, “Ultralow noise miniature external cavity semiconductor laser,” Nature Communications 6, 7371 (2015).
  * [18] B. Li, W. Jin, L. Wu, L. Chang, H. Wang, B. Shen, Z. Yuan, A. Feshali, M. Paniccia, K. J. Vahala, and J. E. Bowers, “Reaching fiber-laser coherence in integrated photonics,” Optics Letters 46, 5201 (2021).
  * [19] D. A. Chermoshentsev, A. E. Shitikov, E. A. Lonshakov, G. V. Grechko, E. A. Sazhina, N. M. Kondratiev, A. V. Masalov, I. A. Bilenko, A. I. Lvovsky, and A. E. Ulanov, “Dual-laser self-injection locking to an integrated microresonator,” Optics Express 30, 17094–17105 (2022).
  * [20] N. G. Pavlov, S. Koptyaev, G. V. Lihachev, A. S. Voloshin, A. S. Gorodnitskiy, M. V. Ryabko, S. V. Polonsky, and M. L. Gorodetsky, “Narrow-linewidth lasing and soliton Kerr microcombs with ordinary laser diodes,” Nature Photonics 12, 694–698 (2018).
  * [21] A. S. Raja, A. S. Voloshin, H. Guo, S. E. Agafonova, J. Liu, A. S. Gorodnitskiy, M. Karpov, N. G. Pavlov, E. Lucas, R. R. Galiev, A. E. Shitikov, J. D. Jost, M. L. Gorodetsky, and T. J. Kippenberg, “Electrically pumped photonic integrated soliton microcomb,” Nature Communications 10, 680 (2019).
  * [22] N. M. Kondratiev, V. E. Lobanov, E. A. Lonshakov, N. Y. Dmitriev, A. S. Voloshin, and I. A. Bilenko, “Numerical study of solitonic pulse generation in the self-injection locking regime at normal and anomalous group velocity dispersion,” Optics Express 28, 38892 (2020).
  * [23] A. S. Voloshin, N. M. Kondratiev, G. V. Lihachev, J. Liu, V. E. Lobanov, N. Y. Dmitriev, W. Weng, T. J. Kippenberg, and I. A. Bilenko, “Dynamics of soliton self-injection locking in optical microresonators,” Nature Communications 12, 235 (2021).
  * [24] G. Lihachev, W. Weng, J. Liu, L. Chang, J. Guo, J. He, R. N. Wang, M. H. Anderson, Y. Liu, J. E. Bowers, and T. J. Kippenberg, “Platicon microcomb generation using laser self-injection locking,” Nature Communications 13, 1771 (2022).
  * [25] W. Jin, Q.-F. Yang, L. Chang, B. Shen, H. Wang, M. A. Leal, L. Wu, M. Gao, A. Feshali, M. Paniccia, K. J. Vahala, and J. E. Bowers, “Hertz-linewidth semiconductor lasers using CMOS-ready ultra-high-Q microresonators,” Nature Photonics 15, 346–353 (2021).
  * [26] W. Liang, D. Eliyahu, V. S. Ilchenko, A. A. Savchenkov, A. B. Matsko, D. Seidel, and L. Maleki, “High spectral purity Kerr frequency comb radio frequency photonic oscillator,” Nature Communications 6, 7957 (2015).
  * [27] G. Lihachev, J. Riemensberger, W. Weng, J. Liu, H. Tian, A. Siddharth, V. Snigirev, V. Shadymov, A. Voloshin, R. N. Wang _et al._ , “Low-noise frequency-agile photonic integrated lasers for coherent ranging,” Nature communications 13, 1–10 (2022).
  * [28] C. Xiang, J. Liu, J. Guo, L. Chang, R. N. Wang, W. Weng, J. Peters, W. Xie, Z. Zhang, J. Riemensberger _et al._ , “Laser soliton microcombs heterogeneously integrated on silicon,” Science 373, 99–103 (2021).
  * [29] J. Guo, C. A. McLemore, C. Xiang, D. Lee, L. Wu, W. Jin, M. Kelleher, N. Jin, D. Mason, L. Chang _et al._ , “Chip-based laser with 1 Hertz integrated linewidth,” arXiv preprint arXiv:2203.16739 (2022).
  * [30] J. Ling, J. Staffa, H. Wang, B. Shen, L. Chang, U. A. Javid, L. Wu, Z. Yuan, R. Lopez-Rios, M. Li _et al._ , “Self-injection-locked second-harmonic integrated source,” arXiv preprint arXiv:2207.03071 (2022).
  * [31] W. Liang, V. Ilchenko, A. Eliyahu, D. Savchenkov, A. Matsko, and D. Seidel, “Ultralow noise miniature external cavity semiconductor laser,” Nat. Commun. 6, 7371 (2015).
  * [32] N. M. Kondratiev, V. E. Lobanov, A. V. Cherenkov, A. S. Voloshin, N. G. Pavlov, S. Koptyaev, and M. L. Gorodetsky, “Self-injection locking of a laser diode to a high-Q WGM microresonator,” Opt. Express 25, 28167–28178 (2017).
  * [33] Z. Zhang, M. Li, J. Guo, B. Du, and R. Zheng, “A portable tunable diode laser absorption spectroscopy system for dissolved CO2 detection using a high-efficiency headspace equilibrator,” Sensors 21, 1723 (2021).
  * [34] P. Laurent, A. Clairon, and C. Breant, “Frequency noise analysis of optically self-locked diode lasers,” IEEE J Quantum Electron. 25, 1131–1142 (1989).
  * [35] M. Gorodetsky, A. Pryamikov, and V. Ilchenko, “Rayleigh scattering in high-Q microspheres,” JOSA B 17, 1051–1057 (2000).
  * [36] A. E. Shitikov, O. V. Benderov, N. M. Kondratiev, V. E. Lobanov, A. V. Rodin, and I. A. Bilenko, “Microresonator and laser parameter definition via self-injection locking,” Physical Review Applied 14, 064047 (2020).
  * [37] M. Gorodetsky and V. Ilchenko, “High-Q optical whispering-gallery microresonators: precession approach for spherical mode analysis and emission patterns with prism couplers,” Optics Communications 113, 133–143 (1994).
  * [38] A. Savchenkov, S. Borri, M. S. de Cumis, A. Matsko, P. D. Natale, and L. Maleki, “Modeling and measuring the quality factor of whispering gallery mode resonators.” Appl. Phys. B 124, 171 (2018).
  * [39] X. Zhang, J. Pouls, and M. C. Wu, “Laser frequency sweep linearization by iterative learning pre-distortion for FMCW LiDAR,” Optics Express 27, 9965–9974 (2019).
  * [40] S. Donati, “Developing self-mixing interferometry for instrumentation and measurements,” Laser & Photonics Reviews 6, 393–417 (2012).
  * [41] L. Lu, J. Yang, L. Zhai, R. Wang, Z. Cao, and B. Yu, “Self-mixing interference measurement system of a fiber ring laser with ultra-narrow linewidth,” Optics Express 20, 8598–8607 (2012).


# Constraining viscous dark energy models with the latest cosmological data

Deng Wang1<EMAIL_ADDRESS>Yang-Jie Yan2<EMAIL_ADDRESS>Xin-He
Meng2<EMAIL_ADDRESS>1Theoretical Physics Division, Chern Institute of
Mathematics, Nankai University, Tianjin 300071, China
2Department of Physics, Nankai University, Tianjin 300071, China

###### Abstract

Based on the assumption that the dark energy possessing bulk viscosity is
homogenously and isotropically permeated in the universe, we propose three new
viscous dark energy (VDE) models to characterize the accelerating universe. By
constraining these three models with the latest cosmological observations, we
find that they just deviate very slightly from the standard cosmological model
and can alleviate effectively the current $H_{0}$ tension between the local
observation by the Hubble Space Telescope and the global measurement by the
Planck Satellite. Interestingly, we conclude that a spatially flat universe in
our VDE model with cosmic curvature is still supported by current data, and
the scale invariant primordial power spectrum is strongly excluded at least at
the $5.5\sigma$ confidence level in three VDE models as the Planck result. We
also give the $95\%$ upper limits of the typical bulk viscosity parameter
$\eta$ in three VDE scenarios.

##  I Introduction

There is no doubt that modern cosmology has entered a new data-driven era.
During the past nearly two decades, a large number of observational evidences
such as Type Ia supernovae (SNIa) 1 ; 2 , baryonic acoustic oscillations (BAO)
3  and cosmic microwave background (CMB) anisotropies 4 ; 5  have confirmed
the late-time acceleration of the universe and verified that our universe
consists approximately of $30\%$ matter including baryons and cold dark matter
(CDM) and $70\%$ dark energy components. The simplest candidate of dark energy
is the so-called cosmological constant scenario, namely the $\Lambda$CDM
model, which has been confirmed as a phenomenologically viable cosmological
model, once again, by the recent Planck CMB measurements with high precision 6
. However, this well-known model suffers from two severe problems 7 : (i) Why
is the value of the cosmological constant unexpectedly small with respect to
any physically meaningful scale, except the current horizon scale ? (ii) Why
this value is not only small, but also surprisingly close to another unrelated
physical quantity, the present-day matter density ? Up to date, due to the
fact that we still know little about the realistic nature of the dark energy,
to resolve these puzzles, a flood of candidates have been proposed by
cosmologists. Generally, they can be divided into two main classes, i.e., the
dark energy models 8 ; 9 ; 10 ; 11 ; 12 ; 13 ; 14 ; 15 ; 16 ; 17 ; 18 ; 19 ;
20 ; 21 ; 22 ; 23 ; 24  and modified theories of gravity (GR) 25 ; 26 ; 27 ;
28 ; 29 ; 30 ; 31 ; 32 . The former case is aimed at introducing a new cosmic
fluid or matter field within the framework of the theory of general
relativity. On the contrary, assuming the breakdown of Einstein’s gravity, the
latter case modifies the standard GR lagrangian based on some concrete
physical mechanism.

Based on the thermodynamics point of view, dissipative process is a universal
property of any realistically physical phenomena. In the literature, as a
consequence, several authors 33 ; 34 ; 35  have naturally used the exotic
fluids from dissipative process, which belongs to the former case mentioned
above, to explain the cosmic acceleration. The concretely physical realization
of dissipative process is the addition of shear and bulk viscosity terms into
the stress-energy tensor of a cosmic fluid 36 ; 37 ; 38 ; 39 ; 40 ; 41 ; 42 ;
43 . Furthermore, assuming the cosmological principle — the isotropic and
homogenous universe, the shear viscosity vanishes and one only needs to
consider the impacts of bulk viscosity on the evolution of the universe in a
given cosmological model. The bulk viscosity effects have been widely
investigated in the cosmological context, for example, inflation in a viscous
fluid model 44 , the gravitino abundance within bulk viscosity cosmology 45 ,
viscous dark matter (VDM) 46 ; 47 ; 48 ; 49 , viscous dark energy (VDE) 14 ;
15 ; 16 ; 17 ; 18 ; 19 ; 50 ; 51 ; 52 ; 53 ; 54 ; 55 , and so forth. Several
sources of viscous fluid can be described as follows: moving cosmic strings
through the comic magnetic fields, magnetic monopoles in monopole interactions
effectively experience various bulk viscosity phenomena 56 ; 57  or primordial
quantum particle productions and their interactions 58 ; 59 ; 60 ; 61 ; 62 ;
63 .

By observing the Einstein equations, one can find that the left-hand side
involves the geometry or gravity, and on the other side the stress-energy
tensor is also related to a quantum field vacuum. Starting directly from the
equation of state (EoS), VDE models are aimed at studying the evolutionary
properties of the universe and explore whether a given EoS of DE can
effectively characterize the observational data. A majority of previous VDE
models are dedicated to single component for the dark sector, i.e., regarding
the DM and DE as a unified dark fluid 14 ; 15 ; 16 ; 17 ; 18 ; 19 ; 23 ; 24 ;
64 . Interestingly, the correspondences between the modified EoS of viscous
cosmological models, the scalar fields, and extended theories of gravity are
also investigated in details 15 ; 65 ; 66 ; 67 . It is noteworthy that the VDE
models have a convenient way to resolve the so-called Big Rip problem 68 ; 69
, which is encountered by the phantom cosmology 13 , since their EoS can
easily exhibit a phantom-barrier-crossing behavior via the bulk viscosity
terms (see v  for a recent review).

In the present work, we focus on proposing three new VDE models and use the
latest cosmological observations including Planck data to constrain them. In
the meanwhile, we also study the abilities of these three models in
alleviating the current Hubble constant $H_{0}$ tension between the indirectly
global measurement by the Planck satellite 70  and the directly local
measurement by the Hubble space telescope 71 . We find that they just deviate
very slightly from the $\Lambda$CDM model and can relieve effectively the
$H_{0}$ tension.

This work is organized as follows. In the next section, we introduce briefly
our three new VDE models. In Section III, we investigate numerically the
effects of bulk viscosity on the evolutional behaviors of the universe in
these models. In Section IV, we describe the latest observational data and our
analysis method. In Section V, we exhibit our results. The discussions and
conclusions are presented in the final section.

##  II VDE models

The dynamics of a Friedmann-Robertson-Walker (FRW) universe are governed by
the following two Friedmann equations

$\frac{\dot{a}^{2}}{a^{2}}=\frac{\rho}{3},$ (1)
$\frac{\ddot{a}}{a}=-\frac{\rho+p}{6},$ (2)

where $a$, $\rho$ and $p$ are the scale factor, energy density and pressure of
the cosmic fluid, respectively, and the dot denotes the derivative with
respect to the cosmic time $t$. Note that we use the units $8\pi G=c=\hbar=1$
throughout this work. Subsequently, we consider a cosmic DE fluid possessing
the bulk viscosity $\zeta$ and its stress-energy tensor is written as

$T_{\mu\nu}^{de}=\rho_{de}U_{\mu}U_{\nu}+\tilde{p}_{de}h_{\mu\nu},$ (3)

where $\rho_{de}$ and $\tilde{p}_{de}$ are the energy density and effective
pressure of VDE fluid, $U_{\mu}=(1,0,0,0)$ the four-velocity of the VDE fluid
in comoving coordinates and $h_{\mu\nu}=g_{\mu\nu}+U_{\mu}U_{\nu}$ the
projection tensor. For a typical VDE fluid, according to Ekart’s theory as a
first order limit of Israel-Stewart scenario with zero relation time, one can
re-express the effective pressure $\tilde{p_{de}}$ at the thermodynamical
equilibrium as $\tilde{p_{de}}=\omega\rho_{de}-\zeta\theta$ [], where
$\omega$, $\zeta$ and $\theta$ is the EoS of perfect DE fluid, bulk viscosity
and expansion scalar, respectively. Furthermore, the corresponding energy
conservation equation for the VDE fluid can be written as

$\dot{\rho_{de}}+\theta(\rho_{de}+\tilde{p_{de}})=0,$ (4)

where $\theta=3H=3\frac{\dot{a}}{a}$ and $H$ is the Hubble parameter that
controls the background evolution of the universe.

We consider the bulk viscosity being proportional to the Hubble parameter,
i.e., $\zeta=\eta H$ in our analysis from first to last. Consequently, the
effective pressure reads $\tilde{p}_{de}=\omega\rho_{de}-3\eta H^{2}$ now. For
the first VDE model, we just consider a one-parameter extension to the
standard six-parameter $\Lambda$CDM cosmology by fixing $\omega=-1$ and only
adding the viscous terms $-3\eta H^{2}$ into the effective pressure. Hereafter
we call it the V$\Lambda$DE model and its effective pressure is
$\tilde{p_{de}}=-\rho_{de}-3\eta H^{2}$. Subsequently, combining Eqs. (1) with
(4), the dimensionless Hubble parameter $E(z)$ of the V$\Lambda$DE model can
be expressed as

$E(z)=\left[\frac{1}{1+\eta}\Omega_{m}(1+z)^{3}+(1-\frac{1}{1+\eta}\Omega_{m})(1+z)^{-3\eta}\right]^{\frac{1}{2}},$
(5)

where $z$ and $\Omega_{m}$ are the redshift and present-day matter density
parameter, respectively. Notice that for the V$\Lambda$DE model, we have
ignored the contribution from the radiation and spatial curvature components.
One can also find that the V$\Lambda$DE model reduces to the standard
$\Lambda$CDM one when $\eta=0$ in Eq. (5). For the second VDE model, we make a
simple one-parameter generalization to the V$\Lambda$DE model by setting the
EoS of DE $\omega$ as a free parameter and we call it the V$\omega$DE model
hereafter. Furthermore, we can exhibit its dimensionless Hubble parameter
$E(z)$ as

$E(z)=\left[\frac{\omega}{\omega-\eta}\Omega_{m}(1+z)^{3}+(1-\frac{\omega}{\omega-\eta}\Omega_{m})(1+z)^{3(1+\omega-\eta)}\right]^{1/2}.$
(6)

For the third VDE model, for the first time, we consider the bulk viscosity
effects on the spatial curvature by adding the curvature contribution into the
cosmic pie. More specifically, we still consider a simple one-parameter
extension to the V$\Lambda$DE model by fixing $\omega=-1$ and we call it the
VKDE model in the following context. The corresponding dimensionless Hubble
parameter $E(z)$ of the VKDE model is

$E(z)=\left[\frac{2}{2+3\eta}\Omega_{k}(1+z)^{2}+\frac{1}{1+\eta}\Omega_{m}(1+z)^{3}+(1-\frac{2}{2+3\eta}\Omega_{k}-\frac{1}{1+\eta}\Omega_{m})(1+z)^{-3\eta}\right]^{1/2},$
(7)

where $\Omega_{k}$ is the present-day curvature density parameter. Note that
the $\Lambda$CDM cosmology recovers in this case when $\eta=\Omega_{k}=0$.
Furthermore, combining the expression of the effective pressure of VDE fluid
$\tilde{p_{de}}=-\rho_{de}-3\eta H^{2}$ and Eq. (4), we obtain, respectively,
the effective EoS of VDE fluid for the V$\Lambda$DE, V$\omega$DE and VKDE
models as

$\omega_{de}(z)=-1-\eta-\frac{\eta\Omega_{m}(1+z)^{3}}{(1-\frac{1-\eta}{1+\eta}\Omega_{m})(1+z)^{-3\eta}-\frac{\eta}{1+\eta}\Omega_{m}(1+z)^{3}},$
(8)
$\omega_{de}(z)=\omega-\eta-\frac{\eta\Omega_{m}(1+z)^{3}}{(1-\frac{\omega+\eta}{\omega-\eta}\Omega_{m})(1+z)^{3(1+\omega-\eta)}+\frac{\eta}{\omega-\eta}\Omega_{m}(1+z)^{3}},$
(9)
$\omega_{de}(z)=-1-\eta-\frac{\eta[\Omega_{m}(1+z)^{3}+\Omega_{k}(1+z)^{2}]}{(1-\frac{2-3\eta}{2+3\eta}\Omega_{k}-\frac{1-\eta}{1+\eta}\Omega_{m})(1+z)^{-3\eta}-\frac{\eta}{1+\eta}\Omega_{m}(1+z)^{3}-\frac{3\eta}{2+3\eta}\Omega_{k}(1+z)^{2}}.$
(10)

Figure 1: The effects of the bulk viscosity coefficient $\eta$ on the
dimensionless Hubble parameter $E(z)$ in the V$\Lambda$DE model. Here we have
assumed $\Omega_{m}=0.3$ and the solid (red), dotted (black), dash-dotted
(orange) and dashed (blue) lines correspond to the cases of $\eta=0$
($\Lambda$CDM), 0.003, 0.03 and 0.3, respectively.

Figure 2: The effects of the perfect EoS of DE $\omega$ on the dimensionless
Hubble parameter $E(z)$ for different values of bulk viscosity coefficient
$\eta$ in the V$\omega$DE model. From left to right, fixing $\Omega_{m}=0.3$,
we consider the cases of $\eta=0.003,0.03,0.3$ and the solid (red), dotted
(black), dash-dotted (orange) and dashed (blue) lines correspond to the case
of $\omega=0$, -1.005, -1.05 and -1.5, respectively.  Figure 3: The effects of
the spatial curvature $\Omega_{k}$ on the dimensionless Hubble parameter
$E(z)$ in the VKDE model. Here we have assumed $\Omega_{m}=0.3$ and
$\eta=0.003$ and the solid (red), dotted (black), dash-dotted (orange) and
dashed (blue) lines correspond to the cases of $\Omega_{k}=0$ ($\Lambda$CDM),
0.003, 0.03 and 0.3, respectively.

Figure 4: Left panel: The effects of the bulk viscosity coefficient $\eta$ on
the effective EoS of VDE $\omega_{de}(z)$ in the V$\Lambda$DE model. The solid
(red), dotted (black), dash-dotted (orange) and dashed (blue) lines correspond
to the cases of $\eta=0$ ($\Lambda$CDM), 0.001, 0.01 and 0.1, respectively.
Middle panel: The effects of the perfect EoS of DE $\omega$ on the effective
EoS of VDE $\omega_{de}(z)$ in the V$\omega$DE model. Fixing $\eta=0.1$, the
solid (red), dotted (black), dash-dotted (orange) and dashed (blue) lines
correspond to the cases of $\omega=-1$, -1.001, -1.01 and -1.1, respectively.
Right panel: The effects of the spatial curvature $\Omega_{k}$ on the
effective EoS of VDE $\omega_{de}(z)$ in the VKDE model. Fixing $\eta=0.1$,
the solid (red), dotted (black), dash-dotted (orange) and dashed (blue) lines
correspond to the cases of $\Omega_{k}=0.001$, 0.005, 0.01 and 0.05,
respectively. Here we have assumed $\Omega_{m}=0.3$.

##  III Qualitative analysis of bulk viscosity effects on the evolution of the
universe

It is very interesting to study the impacts of VDE permeated homogenously at
cosmological scales on the evolution of the universe. To be more specific, we
hope to obtain the qualitative relations between the bulk viscosity
coefficient $\eta$ and two background physical quantities, i.e., the
dimensionless Hubble parameter $E(z)$ and effective EoS of VDE
$\omega_{de}(z)$, in our three VDE models.

From Fig. 1, we find that the larger the bulk viscosity is, the smaller the
dimensionless Hubble parameter is, i.e., the smaller the expansion rate of the
universe is (see also Eq. (5)). However, the effects of bulk viscosity
decrease evidently at low redshifts unlike at high ones, and the evolution of
the universe in the V$\Lambda$DE model behaves very close to that of the
standard cosmology. Note that this conclusion is not only valid for the
V$\Lambda$DE model but also holds true for all the simple VDE models.
Subsequently, we also study the impacts of the perfect freedom $\omega$ on the
expansion rate of the universe in the V$\omega$DE model by fixing
$\eta=0.003$, 0.03 and 0.3, respectively. From Fig. 2, we obtain the
conclusion that, for a universe where $\omega$ is negative enough (e.g.
$\omega=-1.5$) and bulk viscosity is large enough (e.g. $\eta=0.3$), the more
negative the perfect EoS of DE $\omega$ is, the faster the universe expands
(see the right panel of Fig. 2). Furthermore, if taking the spatial curvature
into account, we find that the larger the present-day curvature density
$\Omega_{k}$ is, the faster the universe expands (see Fig. 3). However, this
holds only for a very large spatial curvature. If considering the recent
restriction $|\Omega_{k}|<0.005$ by the Planck Satellite 6 , one can find that
the variation of $\Omega_{k}$ affects hardly the evolution of the universe.
For the VDE models, we are also very interested in the evolutional behaviors
of the effective EoS of VDE. Form Fig. 4, we find that a more negative
effective EoS of VDE $\omega_{de}(z)$ corresponds to a larger bulk viscosity
$\eta$ in the V$\Lambda$DE model, a more negative perfect EoS of DE $\omega$
in the V$\omega$DE model and a larger cosmic curvature $\Omega_{k}$ in the
VKDE model, respectively. Moreover, our three VDE models are very close to the
$\Lambda$CDM one when $z$ approaches zero.

The above-mentioned investigations are just some simple qualitative analysis
about the evolutional behaviors of VDE models. In the next section, we utilize
the latest cosmological observations to constrain these three VDE models and
their corresponding parameter spaces are presented as follows:

$\mathbf{P_{V\Lambda DE}}=\\{\Omega_{b}h^{2},\quad\Omega_{c}h^{2},\quad
100\theta_{MC},\quad\tau,\quad\eta,\quad\mathrm{ln}(10^{10}A_{s}),\quad
n_{s}\\},$ (11) $\mathbf{P_{V\omega
DE}}=\\{\Omega_{b}h^{2},\quad\Omega_{c}h^{2},\quad
100\theta_{MC},\quad\tau,\quad\omega,\quad\eta,\quad\mathrm{ln}(10^{10}A_{s}),\quad
n_{s}\\},$ (12)
$\mathbf{P_{VKDE}}=\\{\Omega_{b}h^{2},\quad\Omega_{c}h^{2},\quad
100\theta_{MC},\quad\tau,\quad\Omega_{k},\quad\eta,\quad\mathrm{ln}(10^{10}A_{s}),\quad
n_{s}\\},$ (13)

where $\Omega_{b}h^{2}$ and $\Omega_{c}h^{2}$ are the present-day baryon and
CDM densities, $\theta_{MC}$ denotes the ratio between the angular diameter
distance and the sound horizon at the redshift of last scattering $z^{\star}$,
$\tau$ is the optical depth due to reionization, $\eta$ is bulk viscosity
coefficient, $\omega$ is the perfect EoS of DE, $\Omega_{k}$ is the present-
day spatial curvature, $\mathrm{ln}(10^{10}A_{s})$ and $n_{s}$ are the
amplitude and spectral index of primordial power spectrum at the pivot scale
$K_{0}=0.05$ Mpc-1. Here $h$ is related to the Hubble constant $H_{0}$ by
$h\equiv H_{0}/(100\mathrm{kms^{-1}Mpc^{-1}})$.

Figure 5: The 1-dimensional posterior distributions on the individual
parameters and $2$-dimensional marginalized contours of the V$\Lambda$DE model
using the data combination CSBCHLW. Figure 6: The 1-dimensional posterior
distributions on the individual parameters and $2$-dimensional marginalized
contours of the V$\Lambda$DE model using the data combination CSBCHLW. Figure
7: The 1-dimensional posterior distributions on the individual parameters and
$2$-dimensional marginalized contours of the V$\omega$DE model using the data
combination CSBCHLW. Figure 8: The $2$-dimensional marginalized contours of
the V$\Lambda$DE (blue), V$\omega$DE (red) and VKDE (green) models in the
planes of $\eta-\tau$, $\eta-n_{s}$, $\eta-H_{0}$, $\eta-\Omega_{m}$, $\eta-
t_{age}$ and $\eta-z_{eq}$ by using the data combination CSBCHLW,
respectively. Figure 9: The 1-dimensional posterior distributions of $\tau$,
$n_{s}$, $H_{0}$, $\Omega_{m}$, $t_{age}$ and $z_{eq}$ in the V$\Lambda$DE
(blue), V$\omega$DE (red) and VKDE (green) models using the data combination
CSBCHLW, respectively.

Figure 10: The $2$-dimensional marginalized contour in the $\eta-\sigma_{8}$
plane and 1-dimensional posterior distribution of $\sigma_{8}$ for the
V$\Lambda$DE (blue), V$\omega$DE (red) and VKDE (green) models using data
combination CSBCHLW.

##  IV Observations and methodology

To study quantitatively the properties of bulk viscosity cosmology, we perform
global constraints on our three VDE models by using the latest cosmological
observations, which are exhibited as follows:

$\bullet$ CMB: The CMB data provides a substantially powerful way of
constraining the background properties for a given cosmological model such as
the matter constituents, topology and late-time evolution of the universe. In
this analysis, we use the Planck 2015 temperature and polarization data 6 ,
which includes the likelihoods of temperature (TT) at
$30\leqslant\ell\leqslant 2500$, the cross-correlation of temperature and
polarization (TE), the polarization (EE) power spectra, and the Planck
low-$\ell$ temperature and polarization likelihood at $2\leqslant\ell\leqslant
29$.

$\bullet$ SNIa: SNIa is a powerful probe of cosmology, especially of EoS of
DE. In 2013, Betoule et al. 72  reported the results of an extensive campaign
to improve the relative photometric calibrations between the SNLS and SDSS SN
surveys, and constructed the largest SNIa sample “ Joint Light-curve Analysis
” (JLA) until now including SNLS, SDSS, HST and several samples of low-$z$ SN.
This sample consists of 740 SNIa data points covering the redshift range
$z\in[0.01,1.3]$.

$\bullet$ BAO: BAO are also excellent geometrical probes and largely
unaffected by errors in the nonlinear evolution of the matter density field
and other systematic uncertainties which may have an impact on other
cosmological data. To break the parameter degeneracies from other
measurements, we employ four BAO data points: the 6dFGS sample at effective
redshift $z_{eff}=0.106$ 73 , the SDSS MGS sample at $z_{eff}=0.15$ 74 , and
the LOWZ at $z_{eff}=0.32$ and CMASS $z_{eff}=0.57$ samples of the SDSS-III
BOSS DR12 sample 75 .

$\bullet$ Cosmic chronometer (CC): The observational Hubble parameter data can
be divided into two classes, i.e., the $H(z)$ data extracted from BAO surveys
and CC data. In general, to obtain $H(z)$ data from BAO observations, one has
to model the redshift space distortions and assume an acoustic scale, both of
which require the assumption of a specific model. As a consequence, this class
of data is model-dependent and cannot be used for constraining a cosmological
model. Nonetheless, the CC data determined by using the most massive and
passively evolving galaxies based on the “ galaxy differential age ” method
are model-independent. We refer the readers to Ref. 76  for detailed
information about CC data (see also Ref. 77 ). Therefore, we use 30 CC data
points in our analysis covering the range $z\in[0.07,1.97]$.

$\bullet$ Lensing (L): Gravitational lensing by the large-scale structure
leaves hints on the CMB temperature and polarization that can be measured in
high angular resolution, low noise observations, such as those from Planck. To
constrain the background parameters better, we also include the Planck lensing
likelihood in our analysis, which has given the most powerful measurement with
a 2.5$\%$ constraint on the amplitude of the lensing potential power spectrum
78 .

Using the above cosmological observations, we adopt the Markov Chain Monte
Carlo (MCMC) method to explore the parameter spaces of our three VDE models.
We modifies the public package CosmoMC 79  and Boltzmann code CAMB 80  to
infer the posterior probability distributions of different parameters. To
perform the Bayesian analysis, we choose uniform priors for different
parameters as shown in Tab. 1. In succession, as studying the effects of
different data combinations on the parameter estimations is beyond the scope
of this study, we shall implement the tightest constraint on these VDE models
using the data combination CMB + SNIa + BAO + CC + L, which is abbreviated as
“ CSBCL ” in the following analysis.

Table 1: The prior ranges and the best-fit values and $1\sigma$ marginalized uncertainties of cosmological parameters in the V$\Lambda$DE, V$\omega$DE, and VKDE models using the data combination CSBCL, respectively. Note that we have quoted the (2$\sigma$) 95$\%$ limits for the parameters which cannot be well constrained. Parameters |  Priors | V$\Lambda$DE | V$\omega$DE | VKDE
---|---|---|---|---
$\Omega_{b}h^{2}$ | $[0.005,0.1]$ | $0.02237^{+0.00017}_{-0.00015}$ | $0.02239^{+0.00017}_{-0.00014}$ | $0.02233^{+0.00015}_{-0.00011}$
$\Omega_{c}h^{2}$ | $[0.001,0.99]$ | $0.11758^{+0.00093}_{-0.00120}$ | $0.1174^{+0.0014}_{-0.0012}$ | $0.1183^{+0.0011}_{-0.0013}$
$100\theta_{MC}$ | $[0.5,10]$ | $1.04099\pm 0.00033$ | $1.04106\pm 0.00027$ | $1.04101\pm 0.00030$
$\tau$ | $[0.01,0.8]$ | $0.0837\pm 0.0037$ | $0.0884\pm 0.0073$ | $0.0807^{+0.0030}_{-0.0034}$
${\rm{ln}}(10^{10}A_{s})$ | $[2,4]$ | $3.0945\pm 0.0071$ | $3.104\pm 0.013$ | $3.0974\pm 0.0075$
$n_{s}$ | $[0.8,1.2]$ | $0.9707^{+0.0048}_{-0.0038}$ | $0.9706\pm 0.0053$ | $0.9695^{+0.0040}_{-0.0029}$
$\eta$ | $[0,1]$ |  $<0.00217$ (2$\sigma$) |  $<0.00375$ (2$\sigma$) |  $<0.00282$ (2$\sigma$)
$\omega$ | $[-2,0]$ | — | $-1.001^{+0.012}_{-0.011}$ | —
$\Omega_{k}$ | $[-0.05,0.05]$ | — | — | $0.0017\pm 0.0018$
$H_{0}$ | $[20,100]$ | $68.24^{+0.61}_{-0.43}$ | $68.34^{+0.55}_{-0.61}$ | $68.81\pm 0.62$
$\Omega_{m}$ |  — | $0.3020^{+0.0054}_{-0.0078}$ | $0.3007^{+0.0080}_{-0.0071}$ | $0.2985\pm 0.0054$
$\sigma_{8}$ |  — | $0.8252^{+0.0039}_{-0.0047}$ | $0.8285^{+0.0045}_{-0.0032}$ | $0.8298\pm 0.0056$
$\chi^{2}_{min}$ |  — | 13715.842 | 13714.371 | 13713.944

##  V Analysis results

Making the use of the tightest constraint we can give, the numerical analysis
results are exhibited in Tab. 1 and Figs. 5-7, which include the best-fit
points and corresponding $1\sigma$ errors (or $2\sigma$ bounds) of individual
parameters and 1-dimensional posterior distributions on the individual
parameters and $2$-dimensional marginalized contours of three VDE models. It
is easy to see that the V$\omega$DE and VKDE models give a similar
$\chi^{2}_{min}$ and they give better cosmological fits than the V$\Lambda$DE
one by the increasements of $\Delta\chi^{2}_{min}=1.471$ and 1.898,
respectively. For three VDE models, we obtain, respectively, the $2\sigma$
upper bounds of the typical parameter $\eta$ 0.00217, 0.00375 and 0.00375 (see
Tab. 1). One can find that the latter two models give larger $2\sigma$ upper
limits than the former one, which maybe ascribed to the enlargement of
parameter spaces in the V$\omega$DE and VKDE models. For the V$\omega$DE
model, we get the perfect EoS of DE $\omega=-1.001^{+0.012}_{-0.011}$, which
is consistent with the prediction $\omega=-1.006\pm 0.045$ of Planck under the
assumption of $\omega$CDM at the $1\sigma$ confidence level (CL) 6 . One can
find that the perfect EoS of DE $\omega$ is still compatible with the standard
cosmology ($\omega=-1$) and the best-fit value still prefers very slightly a
phantom cosmology ($\omega<-1$), even if we consider the bulk viscosity
effects of DE fluid. Furthermore, we also find that the spatial curvature
$\Omega_{k}=0.0017\pm 0.0018$ of our universe in the VKDE model is consistent
with zero at the $1\sigma$ CL, which implies that a spatially flat universe is
supported by current data in the framework of viscous cosmology. Meanwhile,
our result prefers a positive best-fit value corresponding to an open universe
and is in good agreement with the restriction $|\Omega_{k}|<0.005$ of Planck.
Additionally, we conclude that the scale invariant Harrison-Zeldovich-Peebles
(HZP) power spectrum ($n_{s}=1$) 81 ; 82 ; 83  is strongly excluded at the
6.10$\sigma$, 5.55$\sigma$ and 7.63$\sigma$ CL in the V$\Lambda$DE,
V$\omega$DE and VKDE models, respectively, and that their predictive values of
the spectral index $n_{s}$ are well consistent with the Planck evaluation
$n_{s}=0.9655\pm 0.0062$ at the $1\sigma$ CL 6 . Recently, the improved local
measurement $H_{0}=73.24\pm 1.74$ km s-1 Mpc-1 from Riess et al. 2016 71
exhibits a strong tension with the Planck 2015 release $H_{0}=66.93\pm 0.62$
km s-1 Mpc-1 70  at the $3.4\sigma$ CL. Using the data combination CSBCL, we
find that the current $H_{0}$ tension can be effectively alleviated from
3.4$\sigma$ to 2.71$\sigma$, 2.69$\sigma$ and 2.40$\sigma$ in the
V$\Lambda$DE, V$\omega$DE and VKDE models, respectively, and that the
consideration of curvature effect in the viscous cosmology appears to behave a
little better than the other two models.

Utilizing the joint constraint CSBCL, we are also very interested in
investigating the correlations between the typical parameter $\eta$ and other
cosmological parameters in three VDE models. From Fig. 8, we find that the
present expansion rate $H_{0}$ is anti-correlated with the bulk viscosity
coefficient $\eta$, which indicates that the universe expands slower with
increasing viscosity of DE and is well consistent with the above qualitative
analysis (see also Eq. (5) and Fig. 1). Subsequently, we conclude that the
viscosity $\eta$ is positively correlated with the present matter density
$\Omega_{m0}$ and the redshift of radiation-matter equality $z_{eq}$, which
implies that the larger viscous effects of DE, the larger the matter fraction
in the cosmic pie is and the later the time when the radiation and matter
densities equal is. We also find that the optical depth $\tau$ is anti-
correlated with the bulk viscosity $\eta$ in the V$\Lambda$DE and V$\omega$DE
models, which means that the larger the viscosity of DE is, the smaller the
optical depth is. However, this is not the case in the VKDE one, where $\tau$
is positively correlated with $\eta$. Moreover, we find that the age of the
universe increases very slightly with increasing bulk viscosity and the
spectral index appears to decrease very slowly with increasing $\eta$.

To study the details of constrained parameters further and compare
conveniently with each other in three VDE models, we exhibit their
1-dimensional posterior distributions in Fig. 9. We find that the value of the
optical depth of VKDE model is larger than those of the left two models (see
Tab. 1), and the values of the optical depth of these three VDE models are a
little higher than the Planck result $\tau=0.066\pm 0.012$ by using CMB and
several low-$z$ probes 6 . Meanwhile, we conclude that, since the values of
$H_{0}$ and $t_{age}$ of VKDE model has evident shifts with respect to
V$\Lambda$DE and V$\omega$DE ones, the VKDE model gives a larger expansion
rate and a lower age of our universe than the other two ones (see Fig. 8).
Furthermore, three VDE scenarios give similar values of $n_{s}$, $\Omega_{m}$
and $z_{eq}$. In addition, from Fig. 10, we find that the amplitude of the rms
density fluctuations today $\sigma_{8}$ in linear regime is anti-correlated
with the viscosity $\eta$, which indicates that the effect of matter
clustering increases with increasing bulk viscosity of DE. One can also
conclude that the V$\Lambda$DE model gives a lower value of the amplitude of
matter fluctuations than the left two scenarios from the right panel of Fig.
10 (see Tab. 1).

##  VI Discussions and conclusions

Our motivation is to study the effects of bulk viscosity on the evolution of
the universe using the latest cosmological observations. Specifically, we
propose three new VDE models, i.e., V$\Lambda$DE, V$\omega$DE, and VKDE,
investigate their background evolutional behaviors via qualitatively numerical
analysis, and place constraints on them using the data combination CSBCL. For
three VDE models, we obtain the corresponding $2\sigma$ upper limits of the
bulk viscosity coefficient. We find that the V$\omega$DE and VKDE models give
a similar $\chi^{2}_{min}$, while they provide better cosmological fits than
the V$\Lambda$DE one, that these three VDE models can alleviate effectively
the current $H_{0}$ tension from 3.4$\sigma$ to 2.71$\sigma$, 2.69$\sigma$ and
2.40$\sigma$, respectively, and that the consideration of curvature effect in
the VKDE model appears to behave a little better in reliving the $H_{0}$
tension than the other two models. For the V$\omega$DE and VKDE models, we
find that our restrictions on the prefect EoS of DE and the spatial curvature
are well consistent with the predictions by the Planck experiment.

Furthermore, for these three VDE models, we also obtain the following
conclusions: (i) The scale invariance of HZP primordial power spectrum is
strongly excluded, while their constrained values of spectral index are in
good agreement with the Planck result 6 ; (ii) The universe expands slower
with increasing bulk viscosity; (iii) With increasing bulk viscosity, the age
of the universe increases very slightly and the spectral index appears to
decrease very slowly; (iv) The effect of matter clustering increases with
increasing bulk viscosity of DE.

The addition of curvature in the VDE scenario can lead to some interesting
conclusions: (i) The VKDE model can relieve the $H_{0}$ tension a little
better than the other two ones; (ii) The optical depth of VKDE model is anti-
correlated with the bulk viscosity, which is not the case in the left two
scenarios; (ii)The VKDE model gives a larger value of optical depth, a lower
evolutional age and a larger expansion rate of the universe. These conclusions
may imply that there exist deeper correlations between the cosmic curvature
and basic cosmological quantities. Undoubtedly, the spatial curvature is a
direct probe of the spacetime topology, and the detection of a significant
deviation from $\Omega_{k}=0$ would have profound meanings for inflationary
cosmology and fundamental physics. This topic is worth being discussed further
in the future work.

Since our results give very tight restrictions on the bulk viscosity
coefficient $\eta$ in three VDE models, combining with constraints on the
perfect EoS of DE and the spatial curvature, one can easily find that the VDE
models just deviate very slightly from the standard cosmology based on the
current cosmological data (see also Figs. 1-4). We expect that future high-
precision gravitational wave observations can give new information about the
evolution of the universe.

##  VII acknowledgements

We thank B. Ratra and S. D. Odintsov for helpful discussions on cosmology. D.
Wang is grateful to L. Xu, Y. Li, W. Yang, F. Melia and Y. Sun for useful
communications. This work is partly supported by the National Science
Foundation of China.

## References

  * (1) A. G. Riess et al. (Supernova Search Team), Astron. J. 116, 1009 (1998).
  * (2) S. Perlmutter, et al. (Supernova Cosmology Project), Phys. Rev. Lett. 83, 670 (1999).
  * (3) D. H. Weinberg et al., Phys. Rep. 530, 87 (2013).
  * (4) C. L. Bennett et al. (WMAP Collaboration), Astrophys. J. Suppl. Ser. 208, 20 (2013).
  * (5) P. Ade et al. (Planck Collaboration), Astron. Astrophys. 571, A16 (2014).
  * (6) P. Ade et al. (Planck Collaboration), Astron. Astrophys. 594, A13 (2016).
  * (7) S. Weinberg, Rev. Mod. Phy. 61, 1 (1989).
  * (8) R. R. Caldwell, Phys. Lett. B 545 23-29 (2002).
  * (9) Y. Fujii, Phys. Rev. D 26, 2580 (1982).
  * (10) B. Ratra, P. J. E. Peebles, Phys. Rev. D 37, 3406 (1988).
  * (11) S. M. Carroll, Phys. Rev. Lett. 81, 3067 (1998).
  * (12) A. Hebecker, C. Witterich, Phys. Rev. Lett. 86, 3339 (2000).
  * (13) R. R. Caldwell, M. Kamionkovski and N.N. Weinberg, Phys. Rev. Lett. 91, 071301 (2003).
  * (14) X. Meng, J. Ren and M. Hu, Commun. Theor. Phys. 47, 379 (2007).
  * (15) J. Ren, X. Meng, Phys. Lett. B 636, 5 (2006).
  * (16) J. Ren, X. Meng, Phys. Lett. B 633, 1 (2006).
  * (17) M. Hu, X. Meng, Phys. Lett. B 635, 186 (2006).
  * (18) X. Meng, X. Dou, Commun. Theor. Phys. 52, 377 (2009).
  * (19) X. Dou, X. Meng, Adv. Astron. 1155, 829340 (2011).
  * (20) J. A. S. Lima, Phys. Rev. D 54, 2571 (1996).
  * (21) P. Wang, X. Meng, Class. Quant. Grav. 22, 283-294 (2005).
  * (22) A. Kamenshchik, U. Moschella and V. Pasquier, Phys. Lett. B 511, 265 (2001).
  * (23) Q. Zhang, G. Yang, Q. Zou, X. Meng and K. Shen, Eur. Phys. J. C 75, 300 (2015).
  * (24) D. Wang, Y. Yang and X. Meng, Eur. Phys. J. C 77, 263 (2017).
  * (25) S. Nojiri, S. D. Odintsov, Phys. Rev. D 68, 123512 (2003).
  * (26) S. Nojiri, S. D. Odintsov, Phys. Rept. 505, 59-144 (2011).
  * (27) S. Nojiri, S. D. Odintsov, Phys. Lett. B 681, 74-80 (2009).
  * (28) T. Jacobson, Einstein-aether gravity: a status report, PoS QG-PH: 020 (2007).
  * (29) T. Jacobson, Phys. Rev. D 81, 101502 (2010).
  * (30) L. Randall, R. Sundrum, Phys. Rev. Lett. 83, 3370 (1999).
  * (31) L. Randall, R. Sundrum, Phys. Rev. Lett. 83, 4690 (1999).
  * (32) G. R. Davli, G. Gabadadze, M. Porrati, Phys. Lett. B 485, 208 (2000).
  * (33) A. Y. Kamenshchik, U. Moschella and V. Pasquier, Phys. Lett. B 511, 265 (2001).
  * (34) J. C. Fabris, S. V. B. Goncalves and P. E. de Souza, Gen. Rel. Grav. 34, 53 (2002).
  * (35) M. C. Bento, O. Bertolami and A. A. Sen, Phys. Rev. D 66, 043507 (2002).
  * (36) C. Eckart, Phys. Rev. 58, 919 (1940).
  * (37) W. Israel, Annals Phys. 100, 310 (1976).
  * (38) W. Israel, J. M. Stewart, Annals Phys. 118, 341 (1979).
  * (39) M. Szydlowski, O. Hrycyna, Annals Phys. 322, 2745 (2007).
  * (40) A. Avelino, U. Nucamendi, J. Cosmol. Astropart. Phys. 0904, 006 (2009).
  * (41) M. Cataldo, N. Cruz and S. Lepe, Phys. Lett. B 619, 5 (2005).
  * (42) A. Di Prisco, L. Herrera and J. Ibanez, Phys. Rev. D 63, 023501 (2001).
  * (43) I. Brevik and O. Gorbunova, Gen. Rel. Grav. 37, 2039 (2005).
  * (44) K. Bamba, S. D. Odintsov, Eur. Phys. J. C 76, 18 (2016).
  * (45) L. Buoninfante, G. Lambiase, Eur. Phys. J. C 77, 287 (2017).
  * (46) W. Zimdahl, D. J. Schwarz, A. B. Balakin and D. Pavón, Phys. Rev. D 64, 063501 (2001).
  * (47) W. Hiplito-Ricaldi, H. Velten and W. Zimdahl, J. Cosmol. Astropart. Phys. 06, 016 (2009).
  * (48) Phys. Rev. D 82, 063507 (2010).
  * (49) H. Velten, D. J. Schwarz, Phys. Rev. D 86, 083501 (2012).
  * (50) R. Colistete, J. C. Fabris, J. Tossa and W. Zimdahl, Phys. Rev. D 76, 103516 (2007).
  * (51) S. Capozziello, V. F. Cardone, E. Elizalde, S. Nojiri and S. D. Odintsov, Phys. Rev. D 73, 043512 (2006).
  * (52) C. Feng, X. Li, Phys. Lett. B 680, 355 (2009).
  * (53) J. Gagnon, J. Lesgourgues, J. Cosmol. Astropart. Phys. 1109, 026 (2011).
  * (54) A. Sasidharan, T. K. Mathew, Eur. Phys. J. C 75, 348 (2015).
  * (55) G. S. Khadekar, S. Ray and X. Meng, [arXiv: 1611.01023v1].
  * (56) B. L. Hu and D. Pavón, Phys. Lett. B 180, 329 (1986).
  * (57) J. P. Ostriker, A. C. Thompson and E. Witten, Phys. Lett. B 180, 231 (1986).
  * (58) S. Weinberg, Astrophys. J. 168, 175 (1971).
  * (59) J. D. Barrow, Phys. Lett. B 180, 335 (1986).
  * (60) J. D. Barrow, Nucl. Phys. B 310, 743 (1988).
  * (61) V. A. Belinskii, I . M. Khalatnikov, Sov. Phys. JETP 42, 205 (1975).
  * (62) O. Grøn, Astrophys. Space Sci. 173, 191 (1990).
  * (63) R. Maartens, Class. Quant. Grav. 12, 1455 (1995).
  * (64) S. Nojiri, S. D. Odintsov, Phys. Rev. D 72, 023003 (2005).
  * (65) S. Capozziello, S. Nojiri and S. D. Odintsov, Phys. Lett. B 634, 93 (2006).
  * (66) S. Capozziello, S. Nojiri and S. D. Odintsov, Phys. Lett. B 632, 597 (2006).
  * (67) D. Wang, X. Meng, [arXiv: 1604.02951v1].
  * (68) I. Brevik, E. Elizalde, S. Nojiri and S. D. Odintsov, Phys. Rev. D 84, 103508 (2011).
  * (69) P. H. Frampton, K. J. Ludwick and R. J. Scherrer, Phys. Rev. D 84, 063003 (2011).
  * (70) I. Brevik et al., [arXiv: 1706.02543v1].
  * (71) N. Aghanim et al. (Planck Collaboration), Astron. Astrophys. 596, A107 (2016).
  * (72) A. G. Riess et al., Astrophys. J. 826, 56 (2016).
  * (73) M. Betoule et al. (SDSS collaboration), Astron. Astrophys. 568, A22 (2014).
  * (74) F. Beutler et al., Mon. Not. Roy. Astron. Soc. 3017, 416 (2011).
  * (75) A. J. Ross et al., Mon. Not. Roy. Astron. Soc. 835, 449 (2015).
  * (76) A. J. Cuesta et al., Mon. Not. Roy. Astron. Soc. 457, 1770 (2016).
  * (77) M. Moresco et al., J. Cosmol. Astropart. Phys. 05 (2016) 014.
  * (78) D. Wang, X. Meng, Phys. Rev. D 95, 023508 (2017).
  * (79) P. Ade et al. (Planck Collaboration), Astron. Astrophys. 594, A15 (2016).
  * (80) A. Lewis, Phys. Rev. D 87, 103529 (2013).
  * (81) A. Lewis, A. Challinor, and A. Lasenby, Astrophys. J. 538, 473 (2000).
  * (82) E. R. Harrison, Phys. Rev. D 1, 2726 (1970).
  * (83) Y. B. Zeldovich, Mon. Not. Roy. Astron. Soc. 160, 1P (1972).
  * (84) P. J. E. Peebles and J. T. Yu, Astrophys. J. 162, 815 (1970).