\section{Introduction}\label{sec:introduction}

Over the past decade, theoretical studies have established a new paradigm for
the accretion of gas into dark matter halos
to fuel star-formation during
galaxy formation \cite[e.g.,][]{dekel09, keres09, nelson13, vogelsberger13}.
This model predicts that galaxies are
fed by cool `streams' of gas, linked to the surrounding intergalactic medium
(IGM) by a web of cosmic filaments \cite[e.g.,][]{bond96, fukugita98}. These
filaments contain a rich reservoir of nearly pristine gas that drives galaxy formation
and evolution, especially in the early Universe \cite[e.g.,][]{keres05, dekel09,vandevoort11,fumagalli11, dekel13, correa15}.
Nevertheless, this fundamental picture is difficult to test observationally.
Direct imaging of the IGM is crucial for examining this standard paradigm of
galaxy formation and further revealing the IGM-galaxy interactions.

The detection of the IGM in emission was suggested a few
decades ago \cite[e.g.,][]{hogan87,gould96}.
Nevertheless, progress was hindered by the faintness of the IGM emission.
By searching around luminous QSOs, the
expected diffuse emission due to recombination should be
enhanced by a few orders of magnitude
within the densest part of the cosmic web \cite[e.g.,][]{cantalupo05, cantalupo12, kollmeier10}.
 With narrowband
imaging, \citet{cantalupo14, hennawi15} and \citet{cai17} discovered a few sources that are
sufficiently luminous for quantitative analysis of diffuse gas emission.
These sources define the ``Enormous Ly$\alpha$ Nebulae
(ELANe)" \cite[e.g.,][]{cai17b, arrigoni18, cai18, arrigoni19};
they are the extrema of Ly$\alpha$ nebulosities at $z \sim 2-3$, with sizes exceeding the diameters of
massive dark matter halos ($\sim 200$ kpc)  with SB$_{\rm{Ly\alpha}}\ge10^{-17}$ erg s$^{-1}$ cm$^{-2}$
arcsec$^{-2}$ and
Ly$\alpha$ luminosities greater than 10$^{44}$ erg s$^{-1}$.

Recent progress in wide-field
integral-field spectrographs (IFS)
on 8-10m telescopes, including MUSE
and KCWI, provide us an indispensible opportunity to
directly study the IGM/CGM around bright sources
at $z=2-4$ by reaching an unprecedented, low surface brightness (SB) of
 a few $\times 10^{-19}$\,erg s$^{-1}$ cm$^{-2}$ arcsec$^{-2}$. This
makes it possible to study emission from the
circumgalactic medium (CGM) and IGM around bright sources.
{ With VLT/MUSE and using a sample of 17 QSOs}, \citet{borisova16}
reveal that Ly$\alpha$ nebulae with projected sizes exceeding 100 kpc
are ubiquitous for QSOs at $z \gtrsim 3$. \citet{arrigoni19} further confirmed such near-ubiquity
using a sample of 61 QSOs at $z\gtrsim3$.
Nevertheless, a systematic IFS survey has not been performed at $z<3$.
A crucial question is whether one can construct a uniform
sample at a lower redshift (e.g., at $z\sim2$), and probe evolution
in the IGM/CGM across cosmic time.

At $z\approx 2$, \citet{arrigoni16} have conducted deep narrowband
images on 15 $z\approx2.2$ QSOs which have a fainter luminosity
compared with QSO sample at $z\gtrsim3$ \citep{borisova16, arrigoni19}.
 \citet{arrigoni16} did not detect bright nebulae with
 the projected size of $>50$ kpc at Ly$\alpha$ SB of 10$^{-17}$ \sbunit.
 Compared with the ubiquitous Ly$\alpha$ nebulae at $z\approx3$, the narrowband
 results seem to suggest a strong evolution of Ly$\alpha$ emitting cool gas from
 $z=3$ to $z=2$. However, before drawing such a conclusion, we note that the
 following effects may partially yield such a low detection rate
 of Ly$\alpha$ nebulae at $z\approx 2$ (also see \S4.1.3 in Arrigoni Battaia et al. 2019):
(1) The redshift of nebular Ly$\alpha$ may have an offset with the systemic redshift
determined by the Mg {\sc ii}\ emission.
If we use a narrowband with the central wavelength consistent
with the Mg {\sc ii}-determined systemic redshift, then the Ly$\alpha$ emission could reside
outside the narrowband;
(2) The full-width-half-maximum (FWHM) of the
narrowband is still much wider than the FWHM of the nebular Ly$\alpha$ emission.
Thus, the contrast between the point-spread-function (PSF) and the diffuse Ly$\alpha$
emission is higher for the narrowband data comparing to the IFU data, which makes
the PSF subtraction more difficult  for the narrowband data.
{ (3) The surface brightness limit that MUSE reaches is deeper than the
narrowband study which may also yield a higher detection rate of the Ly$\alpha$
nebulae.}

To build a uniform sample for studying IGM emission at $z\approx2$ and directly test
its evolution,
we have conducted a new survey using the Keck Cosmic Web Imager (KCWI)
to search for Ly$\alpha$ nebulae associated with QSOs at $z\approx 2$.
With the KCWI, we can conduct, for the first time, a fair comparison with $z\gtrsim3$ MUSE results,
 by selecting QSOs with similar ionizing luminosities as that of the MUSE samples (Borisova et al. 2016; Arrigoni Battaia et al. 2019).
Our goals are to understand the cool gas budget around
 ultraluminous QSOs at $z\approx2$,
 to study the kinematics throughout QSO halos, and to
 probe the evolution of the CGM around massive halos from
$z=3-2$.

In this paper, we present the first KCWI observations on the
Ly$\alpha$ nebulae associated with ultralumious Type-I QSOs
at $z=2.1-2.3$. We organize this paper as follows. In \S2, we introduce our observations and data reduction.
In \S3, we provide the results of the KCWI observations, perform
optimal extraction of Ly$\alpha$ emission around QSOs,
and study the gas kinematics and morphology.
In \S4, we provide a discussion of the results. The discussion is
based on the observational results from
the KCWI, MUSE, and previous narrowband studies.
Throughout this paper when measuring distances, we refer to physical distances unless
otherwise specified.
We assume a $\rm{\Lambda}$CDM cosmology with
$\Omega_m= 0.3$, $\Omega_{\Lambda}=0.7$ and $h=0.70$.

\section{Observations}\label{sec:observations}

In this section, we provide details on the
KCWI instrument configurations, observations,
data reduction pipeline, and post-processing
after the standard reduction pipeline.

\subsection{KCWI Instrument Configuration}

{Keck Cosmic Web Imager (KCWI) (e.g., Morrissey et al. 2018) is
a general purpose, optical IFS that has been installed on the 10 m Keck II telescope.
KCWI provides seeing-limited imaging from the wavelength range of 3500 -- 5700 \AA,
and the spectral resolution can be configured from $R=1000$ to $R =20000$. The
field of view  is $20''\times33''$ for large slicer, $16''\times20''$ for the medium slicer
and $8''\times20''$ for the small slicer.
KCWI is optimal for a survey of gaseous nebulae
at $z\approx2$ because: (1) KCWI has
a high throughput from $\lambda = 3800 - 5500$\AA, optimal
for probing the Ly$\alpha$, C {\sc iv}, and He {\sc ii}\ lines at $z\sim2$.
(2) KCWI has high spectral-resolution modes ($R>4000$)
which resolve the gas kinematics,
and (3) KCWI has a relatively large
field-of-view (FoV) to  cover extended Ly$\alpha$ nebulae.
Furthermore, KCWI nicely complements the characteristics
of MUSE which thrives at $\lambda > 5000$\AA.  }

Data was taken with the Keck/KCWI instrument between
November 15 2017 and January 30 2019. The seeing
varied in the range of 0.7 -- 1.1 arcsec (FWHM of the
Gaussian at $\approx4000$\AA, measured in the combined
40 min datacubes). The information of the QSO fields
are summarized in Table 1.

For our program, we configured KCWI
with the BM grating and medium slicer which
yields a FoV of 16.8$''$ {perpendicular to slicer}
(24 slicers) and  20$''$ { along the slicer.} {
We also use BM grating and large slicer to observe one of our QSOs: Q1444}

This FoV is sufficient to map the gas around QSO host halos
to a radius of $\approx 100$\,kpc at $z\approx2.3$. {This setting can provide} a spatial sampling
of { $ \frac{20''}{24\ {\rm{slicers}}}\approx 0.67''$} along the slicer and is seeing-limited
perpendicular to the slicer. The spectral resolution
is $R=4000$.  We observed at central wavelengths ranging
from $\lambda= 3900$\AA\ -- 4100\AA\ to cover the Ly$\alpha$ emission of
each QSO in the sample.

The total exposure time for each target is 40 minutes,
which consists of four
10-minute individual exposures. The observing procedure is as follows:
we divide the entire QSO sample into { several sub-groups},
with each sub-group consisting of two to several QSOs
separated by $\lesssim 3$ degrees on the sky.
Also, we require that
the redshift offset $\Delta z$ between each of our QSOs is greater than
$|\Delta z| > 0.05$. This procedure
was taken because the sky is determined using a nearby
offset-target with Ly$\alpha$ emission lines at
different wavelengths (see details in \S2.3.3).

\subsection{KCWI QSO Sample at $z\approx2$}

Our ultraluminous QSO sample at $z\approx2$ is selected from the
SDSS-IV/eBOSS database \cite[e.g.,][]{paris17}, restricted to $2.1<z<2.3$. The lower
limit of $z=2.1$ is set by the blue limit to the sensitivity of KCWI.
The purpose for constraining $z\le 2.3$ is that there
is currently no systematic IFU survey at $z\le2.3$.
Further, at this redshift, H$\alpha$ emission can be observed
from the ground for additional insight into the mechanism(s)
powering Ly$\alpha$ \citep[e.g.][]{leibler18}.
Our KCWI sample can be compared with previous VLT/MUSE observations
at $z>3$ \citep{borisova16,arrigoni19} and Gemini and Keck narrowband
surveys at $z\approx2.3$ \cite[e.g.,][]{arrigoni16}.

We select our QSO sample using the following criteria: (1)
ultraluminous QSOs with $i_{\rm{mag}}<18.5$;
(2) each source must have at
least one other ultraluminous QSOs within $3$\,degree separation,
and with a redshift offset $|\Delta z| >0.05$ (see \S2.3.3 for
details of sky subtraction).
Using these selection criteria, we
select 16 QSOs (listed in Table~\ref{table:QSO_snapshot}).
Note that the SDSS QSO
density peak is at $z\approx2.3$ \cite[e.g.,][]{paris17}, providing us a large database
for selecting targets. The QSOs selected have a median
$i$-band magnitude of 17.7, 0.2 magnitude brighter than the
median $i$-band magnitude of 17 $z\approx3.1$ QSOs in
\citet{borisova16}, and 0.5 magnitude brighter than
the 61 $z\approx3.2$ QSOs
described in \citet{arrigoni19}. We further calculated
$L_{\rm{1450}}$ which is the luminosity $\nu L_{\rm{\nu}}$ at
$\lambda= 1450$\AA. For our KCWI QSO sample, the median
$L_{\rm{1450}}$ is $1.7\times10^{13}$ $L_\odot$. The luminosity of our $z=2$ QSOs,
on average, is similar
to \cite{arrigoni19}.
The QSOs of \cite{borisova16} have a median
luminosity of
$3.2\times10^{13}$ $L_\odot$, i.e.\ $1.9\times$
the median $L_{\rm{1450}}$ of our KCWI QSO sample.
We summarize the QSO luminosity of the
three samples in Figure~\ref{fig:sample}.

\subsection{Data Reduction}
\label{sec:data_reduction}

In this section, we provide a detailed description of the
data reduction, including the standard KCWI pipeline
and our post-processing steps using
the CubeExtractor package
(Cantalupo et al. 2019).

\subsubsection{KCWI Standard Pipeline }

The KCWI pipeline v1.0\footnote{ https://github.com/kcwidev/kderp/}, released in Mar. 17, 2018 was adopted to reduce
our data \citep{morrissey18, cai18}.
For each image, we first subtracted the bias, correcting the pixel-to-pixel
variation using flat-field images,
removing cosmic-rays, and error image creation.
Then, continuum-bar images were used to conduct a
geometric transformation, and the ThAr arc images were analyzed for
wavelength calibration.
At this stage, the datacube was constructed.
Then, the twilight flats were used to correct the
slice-to-slice variance. Each individual image was
  flux calibrated using a spectrophotometric standard star
taken at the beginning of the night.

\subsubsection{CubeFix: Improving the Flat-fielding of the Cubes}

Our scientific goals require the analysis of
diffuse, extended emission at low surface brightness.
Therefore, we must control for systematic
variations across the datacube.
We used custom tools for flat-fielding correction
and sky-subtraction. The procedures
are part of the CubExtractor package \citep{cantalupo19}
which was developed
to improve the detection of faint, low surface brightness emission
in IFU datacubes \citep[e.g. MUSE,][]{borisova16}.
We used the CubeFix routine to correct the systematic
errors due to flat-fielding. Then,
we constructed a medium-band image which is collapsed
using 300 channels in the wavelength direction. The
slice-by-slice correction is then calculated using the medium-band image.
The flat-fielding correction is performed as a self-calibration
using sky-lines and sky continuum as a
uniform source to re-calibrate each individual slice of the IFU.
Sources are masked
in this procedure to minimize the self-calibration errors.
With these steps, the
residual is at a level of less than 0.1\% of the sky.

\subsubsection{Sky Subtraction and Coaddition}

Sky-subtraction was then performed on each individual, flat-field corrected cube
using our custom procedure (CubeSharp routine in the CubExtractor package).
As noted above, we have associated each QSO with another that lies
within three degree separation on the sky.
We estimate the sky, channel-by-channel, from the unsubtracted
datacube of the offset-target.
The sky is calculated by taking the median after masking sources.
In each wavelength, all pixels with 3-$\sigma$ above the median {were} masked as sources, and we repeated this process 10 times to construct a final source mask.
The offset-target always has a redshift offset of $\Delta z > 0.05$ from the main target
({$\Delta z>0.05$ corresponds to $\gtrsim 120$ wavelength channels}).
This criterion insures that
the Ly$\alpha$ emission of the offset-target lies far from the emission from the target so that
we can use the sky determined from the offset-target field around the wavelength of target's Ly$\alpha$ emission.
For each exposure, we determined the QSO center; and then
aligned each exposures according
to this value. We
performed a weighted mean with inverse variance
weighting to construct the final datacube.
The variance images we used are generated from the KCWI pipeline.
{ The seeing condition for each target has been taken into account
in each variance images. The variance amplitude is proportional to the inverse of the
 seeing full-width-half-maximum (FWHM).}

\subsection{Point-spread-function (PSF) and Continuum subtraction}

To search for extended, faint Ly$\alpha$ emission distinct
from the QSO flux,
we applied the following procedure to subtract
the central, QSO emission.
For each wavelength channel, we produce a pseudo-narrowband image
with a width of 300 spectral channels ($\approx150$ \AA).
 For our analysis, we used the median of these $\approx300$
wavelength channels for constructing the pseudo-broadband images. We found that the number of wavelength
channels we adopted provides a good compromise
between capturing wavelength PSF variations and
obtaining a good signal-to-noise ratio in the
pseudo-broadband image.
Then, we rescaled the empirical PSF
according to the integrated flux within the $1''\times1''$
area around the QSO centroid.
Then, we subtracted the PSF at each wavelength
channel \cite[e.g.,][]{herenz15,borisova16}.
We computed the rescaling factor using an averaged-sigma-clip algorithm to
minimize cosmic ray effects.

 We further removed continuum sources in each spaxel (a spectrum in the datacube)
 of the cube using a median-filtering approach to construct
 a continuum image. The window size
 of the median filter is approximately
 300 pixels ($\approx 150$ \AA).   We then subtract the
 continuum image from each wavelength channel to construct a continuum-subtracted
datacube. Sources with flat continua are expected to
be eliminated with this procedure.
 In several cases, stars or background galaxies
 are not completely removed or
 are over-subtracted by this procedure due to the large window size.
  However, these residuals do not affect our results
  because we mask stars
  or galaxies with bright continuum in the cube before extraction.

  In the left panel of Figure~\ref{fig:reduction}, we show the data
  product just after processing by the KCWI standard pipeline.
  In the right panel of Figure~\ref{fig:reduction}, we
  provide the individual data after applying all of
  the post-process reduction described in this section.
  For low surface brightness measurements,
  the sky subtraction is
  crucial. In the final datacube, we checked that of the sky subtraction residual
  at wavelengths outside Ly$\alpha$ emission.
 In our final datacube, the typical 1-$\sigma$ uncertainty is
 $\approx 9.0 \times10^{-19}$ erg s$^{-1}$ cm$^{-1}$ arcsec$^{-2}$ (1\AA\ wavelength bin).
 The root-mean-square (rms) around zero within an aperture
 of $2"$ is $\pm 4.5 \times10^{-19}$ erg s$^{-1}$ cm$^{-1}$ arcsec$^{-2}$
 in the rest-frame wavelength channels between 1255\AA\ -- 1275 \AA,
 a wavelength range that does not contain obvious emission lines from QSOs.
 These results indicate that our sky subtraction and
 PSF subtraction are effective.

\section{Results}

After performing the set of careful reduction steps
described in the last section, the
 1-$\sigma$ flux density uncertainty per voxel\footnote{voxel: a three-dimensional datapoint in an integral field spectrograph datacube.}

 is about 1-$\sigma$ of
$\approx 3.5\times10^{-19}$ erg s$^{-1}$ cm$^{-2}$  \AA$^{-1}$ around the observed wavelength
of 4000\AA. This measurement is consistent with the 1-$\sigma$ error
in surface brightness we reported in the previous section. In \S2.4, we
further checked our sky-subtraction and  PSF subtraction, and we
confirm that the flux density in the sky-subtracted, PSF-subtracted
cubes  is consistent with zero within the rest-frame wavelengths $\lambda=1255$ -- 1275 \AA.
 In this section, we carefully analyze the physical properties of Ly$\alpha$ nebulae around the QSO sample at $z\approx2$.

\subsection{Optimal Extraction of the Ly$\alpha$ Nebulae}
\label{sec:extraction}

We use the three-dimensional, automatic algorithm
 CubExtractor (Cantalupo et al. 2019) to conduct
optimal extraction of extended Ly$\alpha$ emission
in each datacube.
We first smooth the datacubes and
variance using a Gaussian filter with $\approx 1"\times1"$ aperture.
Then, objects are extracted and detected if they contain
at least six connected voxels above a signal-to-noise (S/N)
of two after smoothing both the signal and the variance \footnote{{
The reasons of choosing six voxel limit set is the following:  To compare with previous work
(e.g., Borisova et al. 2016; Arrigoni Battaia et al. 2019), we choose
the number of connected voxel close to the seeing value. The number of pixels that close to
the seeing ($\approx1''$) is $2$ voxels perpendicular to slicer direction (1.3$''$), and $3$ voxels along
the slicer direction (0.9$''$). Thus, we use 6 voxel limit set in the manuscript. Actually, we checked
to use 2 $\times$ 2 voxel smoothing and $3\times4$ voxel smoothing, and we found that the difference of the
surface brightness is less than 1\%. }
}.
The three-dimensional
segmentation masks output by CubExtractor are then used for our analysis.

In Figure~\ref{fig:optimal_extraction}, we present the optimally extracted images of the detected objects in each KCWI
datacube as done by \citet{borisova16}.
Each image has an angular size of 16.8$''\times 20''$ FoV. The images have been
generated by (i) selecting all voxels in the PSF-subtracted and continuum-subtracted KCWI cubes,
(ii) applying the corresponding 3-D masks of each nebula output
by CubExtractor, and then
(iii) integrating the flux along the wavelength direction.
These images are similar to pseudo-NB images, but the width of
the filter is adjusted for each spaxel to include only signal above S/N $>2$.
The width of the pseudo-filter varies from one channel
(typically at the edges of the object) to a few tens of wavelength channels
in the brightest parts of the sources. The white contour in Figure~\ref{fig:optimal_extraction}
represents the 2-$\sigma$ uncertainty of {SB $\approx1.8\times10^{-18}$ erg s$^{-1}$ cm$^{-2}$ arcsec$^{-2}$}.
The SB uncertainty is calculated in 1 arcsec$^2$ area and in a wavelength bin of 1\AA\ (also see \citet{arrigoni19}).
Note the 2-$\sigma$ SB limit is slightly different field-by-field.
We summarize the SB limit for each source
in Table~1.

\subsection{Detection rate of giant Ly$\alpha$ nebulae}

Above our KCWI detection limit of
$\approx 1.8\times 10^{-18}$ erg s$^{-1}$ cm$^{-1}$ arcsec$^{-2}$,
we find that {14 out of 16 QSOs} in our sample
have a detected nebula with diameter $\gtrsim50$ kpc.
Among them, Q1227, Q1228, Q1230, and Q1416 are
enormous Ly$\alpha$ nebulae at $z\approx2$ with
the SB$_{\rm{Ly\alpha}}>10^{-17}$ erg s$^{-1}$ cm$^{-2}$ arcsec$^{-2}$
for $\ge100$ kpc,
and their projected sizes exceed the FoV of KCWI.
Q0048 and Q1426 fields contain compact Ly$\alpha$ emitters in the
KCWI fields, possibly indicating a strong overdense nature in these fields.

\subsection{Surface brightness of the Ly$\alpha$ Emission}

In this section, we measure and
present the radial Ly$\alpha$ profiles of the nebulae.
For comparison with previous works,
we use circular-averaged
surface brightness (SB) profiles centered on the QSO continua.
The SB profile is calculated using the pseudo-narrowband images.
we integrate over a fixed velocity range
of $\pm 1000$ km s$^{-1}$ around the centroid of Ly$\alpha$
nebular emission to calculate the surface brightness. This allows
us to properly compare the SB calculated in our sample with
previous work \citep[e.g.,][]{borisova16, arrigoni19}.

The individual,
circularly-averaged SB profiles for each QSOs are shown
as light gray lines in Figure~\ref{fig:SB_r75} and Figure~\ref{fig:SB_r90}.
Further, we calculated the median SB using
the full KCWI QSO sample at $z\approx2.2$ (shown in the thick red color
in the Figure~\ref{fig:SB_r75} and Figure~\ref{fig:SB_r90}).
{ In Figure~\ref{fig:annuli}, we show the annulus (white circles)
which are used to calculate the surface brightness in Figure~\ref{fig:SB_r75}
and Figure~\ref{fig:SB_r90} for each QSO.
We marked the centers of each QSO position using the filled
black circle. }
We found that the median Ly$\alpha$ SB can be described
by the following power-law profile centered at the QSO and
valid on the radius of $r\approx15-70$ kpc:

 \begin{equation}
 \begin{split}
 SB_{\rm{KCWI}}(z\approx2.3) = & 3.7 \times 10^{-17} \times (r/40\ \rm{kpc})^{-1.8} \\
 &  \rm{erg} \ \rm{s^{-1}}   \ \rm{cm^{-2}}  \ \rm{arcsec^{-2}}
  \end{split}
 \end{equation}

 Borisova et al. (2016) have conducted a MUSE snapshot survey
 on 17 QSOs at $z\ge 3.1$,
 and found that all of them are associated with
 large Ly$\alpha$ nebula on a spatial extent of $\ge 100$ kpc.
 Here, we denote the SB of Borisova et al. (2016) as SB$_{\rm{B}}$.
 The median SB$_{\rm{B}}({\rm{Ly\alpha}})$
 can be described by the following equation:
 \begin{equation}
 \begin{split}
 SB_{\rm{B}}(z\approx3.1) = & 3.2 \times 10^{-17} \times (r/40\ \rm{kpc})^{-1.8} \\
 &  \rm{erg} \ \rm{s^{-1}}   \ \rm{cm^{-2}}  \ \rm{arcsec^{-2}}
  \end{split}
 \end{equation}

Using a larger sample of 61 QSOs at $z\gtrsim3.2$, \citet{arrigoni19} also found that
the stacked QSO profiles can be fitted with an exponential profile with
$SB(r) = C_{\rm{e}} {\rm{exp}}(-r/r_{\rm{h}})$
or a power law $SB(r) = C_{\rm{p}} r^\alpha$.
These QSOs are characterized by a median redshift of  $z= 3.17$ ($3.03 < z < 3.46$), absolute
$i$ magnitude in the range $-29.67 \le M_i \le  -27.03$ similar
as that of Borisova et al. (2016).
Here, we denote $SB_{\rm{A}}$ as
the median surface brightness derived from \citet{arrigoni19} sample.
 \citet{arrigoni19} found that the
obtained stacked profiles may be better fit by an exponential profile with scale length of
$r_{\rm{h}} = 15.2\pm0.5$ kpc for radio-quiet objects, i.e.,
\begin{equation}
 \begin{split}
SB_{\rm{A}}(z\approx3.1)= & 5.4\times10^{-17} \ \rm{exp(-r/15.2 {\rm{kpc}}}) \\
&  \rm{erg} \ \rm{s^{-1}}   \ \rm{cm^{-2}}  \ \rm{arcsec^{-2}}
  \end{split}
\end{equation}

Here, if we correct for cosmological SB dimming of MUSE \citep{borisova16, arrigoni19}
from $z\approx 3.1$ to $z=2.2$, then the redshift-scaled median SB of the MUSE observations are:

 \begin{equation}
 \begin{split}
 SB^{\rm{scaled}}_{\rm{B}} (z\approx2.3)\approx & 7.6 \times10^{-17} \times (r/40\ \rm{kpc})^{-1.8} \\
  &  \ \rm{erg}   \ \rm{s^{-1}}   \ \rm{cm^{-2}}  \ \rm{arcsec^{-2}}
\end{split}
 \end{equation}

\begin{equation}
 \begin{split}
SB^{\rm{scaled}}_{\rm{A}}(z\approx2.3)\approx & 13.7\times10^{-17} \ \rm{exp(-r/15.2 {\rm{kpc}}}) \\
&  \rm{erg} \ \rm{s^{-1}}   \ \rm{cm^{-2}}  \ \rm{arcsec^{-2}}
  \end{split}
\end{equation}
where the $SB^{\rm{scaled}}$ means that we corrected the surface brightness
profile of both \citet{borisova16} ($SB_{\rm{B}}$)
and \citet{arrigoni19} ($SB_{\rm{A}}$) from $z\approx3.1$
to $z\approx2.3$, scaled by the cosmological surface brightness dimming of the $(1+z)^4$ factor.
{ We present the scaled data as the purple and orange points }
in Figure~\ref{fig:SB_r75} and Figure~\ref{fig:SB_r90}.
Note the error bar is not just the statistical error, and the error bar indicates the 25\% -- 75\% percentile
 in
the Figure~\ref{fig:SB_r75}; and 10\% -- 90\% percentile in the Figure~\ref{fig:SB_r90}
around the median surface brightness of $z\approx3$ QSO samples.
The median SB of \citet{arrigoni19} is consistent with \citet{borisova16}.
Note the median re-scaled SB of \citet{arrigoni19} is
a factor of $1.5\times$ that of \citet{borisova16} over $20 - 40$
kpc in radius, and one can also see this using Eqs.(4) and (5).
Both MUSE samples suggest that the re-scaled median SB profiles at $z\approx3$ is
a factor of $2-3\times$ our KCWI results at $z\approx2$.

\begin{figure}[h]
\includegraphics[width=0.5\textwidth,height=0.28\textheight]{L_1450_distribution.pdf}
 \caption{The distribution of the QSO luminosity at the
 rest-frame of 1450 \AA\ at $z\approx2$ (Red line).
 The $z\approx2$ QSO sample has a luminosity of  $L_{\nu, \rm{1450\AA}}= 10^{31.48\pm0.32}$ erg s$^{-1}$ Hz$^{-1}$.
 The blue line represent the QSOs at $z\approx3$ from
 Borisova et al. (2016) which { have} a QSO luminosity of
 $L_{\nu,\rm{1450\AA}}= 10^{31.77\pm0.25}$ erg s$^{-1}$ Hz$^{-1}$.
 The yellow dot-dashed line represents
 the QSO luminosity at 1450 \AA\ at $z\approx3$ in
 Arrigoni Battaia et al. (2019) which has a QSO luminosity of
  $L_{\nu,\rm{1450\AA}}= 10^{31.49\pm0.23}$ erg s$^{-1}$ Hz$^{-1}$.
 Our KCWI QSO sample has a similar {median bolometric} luminosity as
 Arrigoni Battaia et al. (2019) and is $0.29$ dex fainter
 than Borisova et al. (2016). }
 \label{fig:sample}
\end{figure}

\begin{figure}[h]
\includegraphics[width=0.5\textwidth,height=0.29\textheight]{data_reduction.pdf}
 \caption{Left panel shows a white-light image obtained from a cube reduced with the KCWI standard pipeline
 (v1.0).  Middle cube indicates
 the data product after running our post-process script CubeExtractor (CubeFix routine).
 Note that the scattered light is reduced.
 Right cube shows the Point-spread-function (PSF) and continuum-subtracted
 cube. The images shown here are constructed  using the median of 2000 km s$^{-1}$
 velocity channels between the Ly$\alpha$ and NV region, where we expect no significant line
 emission in this velocity range, and the median flux is consistent with zero. }
  \label{fig:reduction}
\end{figure}

In Figure~\ref{fig:SB_r75} and Figure~\ref{fig:SB_r90}, we also include Arrigoni
Battaia et al. (2016) results, { shown as the dashed
orange curve with error bars}. Using narrowband observations on
QSOs at $z\approx2$, \citet{arrigoni16} obtained the median Ly$\alpha$
SB at $z=2.2$, similar to the redshift of KCWI QSO sample. The QSOs used
in \citet{arrigoni16} are 1.15 magnitude fainter in $i$-magnitude
than our KCWI QSO sample. Interestingly, the narrowband Ly$\alpha$ SB is one
order of magnitude lower than the KCWI observations.
Note the narrow-band sample \citet{arrigoni16} indeed detected extended Ly$\alpha$
in 47\% of the fields with maximum projected sizes of $\lesssim50$ kpc
above $10^{-17}$ erg s$^{-1}$ arcsec$^{-2}$. These fields are consistent with
seven QSOs in our KCWI sample. Nevertheless, our KCWI {survey} detected a larger
fraction of bright nebulae compared to the narrowband
probes, yielding a much higher median SB profile.
In \S4, we will discuss several possible reasons for this discrepancy.
In Figure~\ref{fig:SB_r75} and Figure~\ref{fig:SB_r90}, we also plot the Ly$\alpha$
profile of Lyman break galaxies (LBGs)
at $z\approx3$ \cite[e.g.,][]{steidel11} and Ly$\alpha$ emitters (LAEs) \cite[e.g.,][]{wisotzki16}.
From the comparison, the Ly$\alpha$ profile powered by galaxies are
much fainter than that powered by QSOs.

In Figure~\ref{fig:2Dspec}, we show the two-dimensional (2D) and
the corresponding one-dimensional spectra using
a pseudo-slit  for each QSO. The pseudo-slit for each nebula
is the white contour shown in Figure~\ref{fig:optimal_extraction}.
Each corresponding 1D spectrum is obtained by integrating
all spatial pixels within the white contour in Figure~\ref{fig:optimal_extraction}.
In each sub-figure, we detect high S/N Ly$\alpha$ emission line.
The extended Ly$\alpha$ emission for each QSO have
 FWHM of $400 - 800$ km s$^{-1}$, i.e.\ much narrower than
the QSO Ly$\alpha$ emission (blue spectra in Figure~\ref{fig:2Dspec}).
This confirms that the extended
nebular emission is unrelated to PSF
subtraction residuals.

\begin{figure*}[h!]
\includegraphics[width=1.0\textwidth,height=0.74\textheight]{optimal_extraction.pdf}
 \caption{``Optimally-extracted" Ly$\alpha$ images from PSF and continuum-subtracted KCWI datacubes for each QSO.
 The white bar indicates a physical scale of
 100 kpc. Each image has a size of $16''$ in x-axis and $20''$ in y-axis.
  The images have been produced by collapsing the datacube voxels associated with
  the CubExtractor
 3-D segmentation maps (the ``3D-mask") along the wavelength direction (see \S2).
 The SB is calculated using the CubExtractor 3D segmentation mask
 that includes a different number of spectral resolution bins.
 The 3D segmentation masks
 have been obtained with a signal-to-noise ratio (SNR) threshold of 2 per smoothed voxel,
 similar to \citet{borisova16} and \citet{arrigoni19}.

 {The white thick contours in each image
 corresponds to 2-$\sigma$, 5-$\sigma$ and 10-$\sigma$ SB, with 1-$\sigma_{\rm{SB}} = 5\times10^{-19}$ erg s$^{-1}$
 cm$^{-2}$ arcsec$^{-2}$}. The black dots represent the positions of the central QSOs.}
  \label{fig:optimal_extraction}
\end{figure*}

\begin{figure*}[h!]
\includegraphics[width=0.9\textwidth,height=0.6\textheight]{SB_r_75percent.pdf}
 \caption{Ly$\alpha$ surface brightness (SB) profiles (circularly averaged) as a function of
 radius around ultraluminous QSOs. All errorbars
 represent the 25 -- 75 percentile
 of Ly$\alpha$ SB in each QSO sample. The gray line with
 cyan points represent the SB profile
 for individual QSO at $z\approx2$. Same as \citet{borisova16} and \citet{arrigoni19},
 the SB is calculated in the pseudo-narrowband images with the width of 2000 km s$^{-1}$, centered on
 the nebular Ly$\alpha$ emission (\S3.3).
 The thick red
 curve represents the median Ly$\alpha$
 profile of the KCWI sample at $z\approx2$. The solid line
 is the median of the actual data, and the dashed red
 curve indicates the extrapolation results.
 The red region represents the SB within the
 25\% and 75\% percentile of our KCWI QSO sample.
 Purple represents the median Ly$\alpha$ SB
 profile of 17 QSOs
 at $z\approx3$ (Borisova et al. 2016)
 {while the error bars represent} the 25 and 75 of
 percentile of this sample.
 The solid orange with error bar represent
 the median and 25 -- 75 percentile Ly$\alpha$ SB of 61 QSOs
 at $z\approx3.2$ reported by Arrigoni Battaia et al. (2019).
 Orange points with dashed error bar
 represent the SB of Ly$\alpha$ emission using
 a narrowband filter (Arrigoni Battaia et al. 2016).
 Cyan represents the Ly$\alpha$
 SB for Lyman break galaxies (LBGs)
 at $z\approx3$ (Steidel et al. 2011),
 and blue shows the Ly$\alpha$ profile
 for $z\approx3$ Lyman alpha emitters (LAE)
 (e.g., Wisotzki et al. 2016). }
 \label{fig:SB_r75}
\end{figure*}

\begin{figure*}[h]
\includegraphics[width=0.9\textwidth,height=0.6\textheight]{SB_r_90percent.pdf}
 \caption{Similar to Fig.~\ref{fig:SB_r75}, but showing 10 -- 90 percentile of
 the Ly$\alpha$ surface brightness of each QSO sample.}
 \label{fig:SB_r90}
\end{figure*}

\begin{figure*}[h!]
\includegraphics[width=0.98\textwidth,height=0.7\textheight]{optimal_extraction_annuli.pdf}
 \caption{In this figure, we present the annuli (white circles)
which are used to calculate the surface brightness in Figure~\ref{fig:SB_r75}
and Figure~\ref{fig:SB_r90} for each QSO. Also, the black
filled circle in the center represents the QSO position which is { used to calculate the circularly
averaged SB profiles.}  }
 \label{fig:annuli}
\end{figure*}

\begin{figure*}[h!]
\includegraphics[width=1.0\textwidth,height=0.95\textheight]{2D_spectra.pdf}
 \caption{The 2-D and 1-D spectra of Ly$\alpha$ nebulae at $z\approx2$
 obtained from KCWI datacubes for each ultraluminous QSOs.
 The upper panel shows the 2-D spectra, extracted using the pseudo-aperture defined by the thick white contours defined in Figure~\ref{fig:optimal_extraction} and integrating along the spatial x-axis direction.
 The lower panels show the 1-D spectra of the nebulae (black lines) obtained by integrating the 1-D
 spectra along the spatial direction between two horizontal orange lines in the upper panel.
 We overlay the one-dimensional spectrum of the QSO. Most of radio-quiet nebulae show a Ly$\alpha$ spectral shape (black) very different from that of the QSO (blue),
 confirming that the detected emission is not an artifact of the QSO PSF subtraction. {We overplot as a vertical purple line the expected Ly$\alpha$ velocity from the systemic redshift of the QSO as determined by the SDSS pipeline fitting routine (\citealt{paris17}, also see Table~1). The orange vertical line shows instead the systemic redshift obtained by a fit of the quasar broad Ly$\alpha$ line}.}
 \label{fig:2Dspec}
\end{figure*}

\begin{figure*}[h!]
\includegraphics[width=0.95\textwidth,height=0.5\textheight]{NB_zsys_vs_zLya1.pdf}
 \caption{Comparison between two pseudo-narrowband
  images of Q1228 {obtained with different filter central wavelengths}.
  { The left panel shows a pseudo-narrow band image obtained assuming as the
  filter central wavelength the systemic redshift determined using the MgII emission (Table~1)}.
  {The right image uses instead the nebula Ly$\alpha$ redshift as the filter central wavelength.}
  The narrowband we use has a similar profile
  of NB3950 (e.g., Cantalupo et al. 2014), and the surface brightness is calculated by integrating
  the flux over a wavelength bin of 30 \AA.  }
   \label{fig:compare}
\end{figure*}

\begin{figure*}[h!]
\includegraphics[width=0.7\textwidth,height=0.44\textheight]{morphology.pdf}
 \caption{Plot of the asymmetry $\alpha$,
 i.e. the ratio between the semiminor
 axis b and semimajor axis a,
 versus the area enclosed by the 2$\sigma$ isophote.
 Black indicates
 the QSOs in Arrigoni Battaia et al. (2019) at $z\approx3.1$; and the red represents QSOs in our sample at $z\approx2.2$. Our sample has a median $\alpha$ value of 0.54 (red horizontal dotted line), comparing with the median $\alpha$ of $0.71$ at $z\approx3$ (black horizontal dotted line). Although a larger sample is needed, our current data tentatively suggests that the QSOs at $z\approx2$ could be more asymmetric than that at $z\approx3$.
 }
 \label{fig:alpha_size}
\end{figure*}

\begin{figure*}[h!]
\includegraphics[width=1.0\textwidth,height=0.7\textheight]{optimal_velocity.pdf}
 \caption{This figure shows the flux-weighted velocity map
 for individual QSOs at $z\approx2$ using KCWI. }
 \label{fig:velocity}
\end{figure*}

\begin{figure*}[h!]
\includegraphics[width=1.0\textwidth,height=0.7\textheight]{optimal_velocityD.pdf}
 \caption{This figure shows the flux-weighted velocity dispersion map
 for individual QSOs at $z\approx2$ using KCWI. }
 \label{fig:velocityD}
\end{figure*}

\begin{figure*}[h]
\includegraphics[width=0.9\textwidth,height=0.25\textheight]{L_Lya_FWHM.pdf}
 \caption{The left panel shows that the FWHM of Ly$\alpha$ emission as a function of nebular Ly$\alpha$ luminosity.
 The right panel indicates the size of the FWHM
 of Ly$\alpha$ emission as a function of ionization luminosity. Due to the extremely high
 SNR of Ly$\alpha$ emission, the luminosity measurement has a typical error of $\le1\%$
 and the FWHM measurement has a typical error of $10\%$ (see Table~1).
 No strong
 correlation between the FWHM and the nebular
 Ly$\alpha$ luminosity { is found}. }
 \label{fig:LLya_FWHM}
\end{figure*}

\begin{figure}[h!]
\includegraphics[width=0.5\textwidth,height=0.29\textheight]{LLya_Lnu.eps}
 \caption{The relation between the Ly$\alpha$ luminosity and the QSO ionizing flux $L_{\rm{\nu}LL}$. From this figure, we can see that there is no
 obvious correlation between $L_{\rm{Ly\alpha}}$ and  $L_{\rm{\nu}LL}$, suggesting the nebular powering
 mechism cannot be dominated by optically thick photoionization scenario. The typical statistical
 error of the data points in this figure is less than 1\% (see Table~1).}
 \label{fig:LLya_Lnu}
\end{figure}

\subsection{Morphology of the nebulae}
\label{sec:morphology}

From Figure~\ref{fig:optimal_extraction}, each nebula has a different morphology and size. All of them have asymmetric morphology.
We use the method of \citet{arrigoni19}
to quantify the morphology.
Our results can be used to quantify the morphological  evolution of the Ly$\alpha$ nebulae from $z=3$ to $z=2$.
We used the asymmetry of $\alpha$ to
quantify the morphology, and $\alpha$ can be calculated from the following formulae:

\begin{equation}
\begin{split}
M_{\rm{xx}}\equiv \left\langle\frac{(x-x_{\rm{Neb}})^2}{r^2} \right\rangle_f; \ \
M_{\rm{yy}}\equiv \left\langle\frac{(y-y_{\rm{Neb}})^2}{r^2}\right\rangle_f;  \\
M_{\rm{xy}}\equiv \left\langle\frac{(x-x_{\rm{Neb}})(y-y_{\rm{Neb}})}{r^2}\right\rangle_f
\end{split}
\end{equation}
where $x_{\rm{Neb}}$, $y_{\rm{Neb}}$
are the flux-weighted centroid
for each nebulae within the 2-$\sigma$ isophoto in the 2D
image shown in Figure~\ref{fig:optimal_extraction}, and $r$ is the distance of a point (x,y) from the flux-weighted centroid. The subscript $f$
represents the flux-weighted centroid.
\begin{equation}
Q\equiv M_{\rm{xx}} - M_{\rm{yy}}, \ \ U\equiv 2 M_{\rm{xy}}
\end{equation}
where M$_{\rm{xx}}$, M$_{\rm{xy}}$, M$_{\rm{yy}}$ are called
second-order moments (``Stokes parameters") (see Arrigoni Battaia et al. 2019).
Using Eq.~(6) and Eq.~(7), we define the asymmetry $\alpha$ using the following equation:

\begin{equation}
\alpha = b/a = \frac{(1- \sqrt{Q^2+U^2})}{1+ \sqrt{Q^2+U^2}}
\end{equation}
We find that at $z \approx 2 $, the median of the asymmetry
of our 16 QSOs is 0.54, with a 1$\sigma$ scatter of 0.18 (see Figure~\ref{fig:alpha_size}). Note the
surface brightness for yielding this results are 2$\sigma$ limit of $\approx 2\times10^{-18}$ erg s$^{-1}$ cm$^{-2}$
arcsec$^{-2}$.
The asymmetry of the nebular morphology does not have a strong dependence of the size of the nebulae.
The bright nebulae with
diameters $>100$ kpc (e.g., Q2123, Q1230, Q1416) show clumpy and filamentary structures.
Q2121 shows two major components. Q0048 shows a strongly asymmetric Ly$\alpha$ distribution and
a compact Ly$\alpha$ emitter in the field.
For a comparison, \citet{borisova16}
suggest that most nebulae
at $z>3$ with modest Ly$\alpha$ emission have more
symmetric and circular emission; nebulae with scales
$>200$ kpc show evidence of filamentary structure
and complex multiple components. \citet{arrigoni19} also found that
most of the nebulae at $z>3$ have symmetric and round morphologies,
with a median asymmetry of $\alpha=0.71$. One of the most asymmetric
structure in the sample is the ELAN in \citet{arrigoni19} with $\alpha\approx0.5$.
The $z\approx2$ QSOs have a smaller median asymmetry. Although larger sample is required,
$z\approx2$ QSOs may be more asymmetric and clumpy
compared to the $z\approx3$ samples \citep{arrigoni19}. We will discuss
the implications of these observations in \S4.5.

\subsection{Kinematics of the Ly$\alpha$ Nebula}

Although the Ly$\alpha$ line may be broadened by radiative transfer effects,
 the relative comparison between the kinematics of
different objects is still informative.
In Figure~\ref{fig:velocity}, we show the map of the first moment (flux-weighted velocity)
for each Ly$\alpha$ nebula, centered on
the peak of the integrated Ly$\alpha$ emission of each nebula.
Q1416, Q1228, may show possible rotation kinematics in
 disk-like structures, possibly
suggesting kinematics of gas inflow predicted by simulations
\cite[e.g.,][]{stewart13}.
Nevertheless, the majority of the nebulae do not show any clear evidences
 of rotation or well-regulated kinematic patterns.

In Figure~\ref{fig:velocityD}, we present the velocity dispersion map
(the second moment of the flux
distribution).
The median value of the velocity dispersion
for each sources ranges from 83 -- 381 km s$^{-1}$ \footnote{Note the FWHM of
the emission line if Gaussian
is approximately a factor of 2.35 $\times$
the velocity dispersion. }.
Our current data does not suggest
any clear correlation between the luminosity and
the velocity dispersion of Ly$\alpha$ emission (Figure~\ref{fig:LLya_FWHM}).
The velocity dispersion shown in Figure~\ref{fig:velocityD} is consistent
with the FWHM calculated from the integrated spectra shown in the Figure~\ref{fig:2Dspec}.
The velocity dispersion of $z\approx2$ QSOs have a consistent range
with that of $z\approx3$ \citep{borisova16, arrigoni19}.
The motions within Ly$\alpha$ nebulosities have amplitudes
consistent with gravitational motions expected in dark matter halos hosting
QSOs ($M_{\rm{halo}}\sim10^{12.5}\ M_\odot$) at these redshifts \cite[e.g.,][]{arrigoni19}.

\section{Discussions}
\label{sec:dicussion}

From the analysis presented in the \S3, it is
clear that our KCWI observations
reveal extended Ly$\alpha$ emission on the projected
scales exceeding 50 kpc around 14 out of 16 of
the $z\approx2$ QSOs (except Q0848 and Q2125).
Further, the median surface brightness
of Ly$\alpha$ in the KCWI study is higher than that of
previous narrowband  studies at similar redshift.
 KCWI results suggest that the typical
circularly-averaged SB profile of Ly$\alpha$ nebulae
at $z\approx2$ is a factor of $2-3$ lower than that at $z\approx3$ determined
by MUSE observations \cite[e.g.,][]{borisova16,arrigoni19}.
In this section, we discuss the implications of the observations,
relying only on the information
coming from the KCWI observations.

 \subsection{A much higher detection rate of Ly$\alpha$ nebulae at $z\approx2$
 compared to narrowband surveys}

The results of the KCWI survey at $z\approx2$ is in contrast with previous
narrowband studies at the same redshift.
In our KCWI survey,
{ 14 out of 16 QSOs} at $z\approx2$ are
associated with Ly$\alpha$ nebulae having
projected linear scales of $>50$ kpc,
with the 2-$\sigma$ SB limit of
$2\times10^{-18}$ erg s$^{-1}$ cm$^{-2}$ arcsec$^{-2}$.
The median profile is
$SB_{\rm{Ly\alpha}}= 3.7 \times 10^{-17} \times (r/40\ \rm{kpc})^{-1.8}$
erg s$^{-1}$ cm$^{-2}$ arcsec$^{-2}$ (see \S3).

The KCWI survey reaches a typical 2-$\sigma$ SB of $1.6\times10^{-18}$ erg s$^{-1}$
cm$^{-2}$ arcsec$^{-2}$ (1 \AA\ wavelength bin),
deeper than previous narrowband surveys (typically
$\approx4\times10^{-18}$ erg s$^{-1}$ cm$^{-2}$ in \citet{arrigoni16}).
The luminous, large Ly$\alpha$ emission
detected in our KCWI survey
is rarely found in previous
narrowband survey
We use dashed yellow error bars in (Figure~\ref{fig:SB_r75} and Figure~\ref{fig:SB_r90})
to represent a large narrowband survey
(Arrigoni Battaia et al. 2016) which study
suggests that the Ly$\alpha$ SB level is $5\times10^{-19}$ erg s$^{-1}$ cm$^{-2}$ arcsec$^{-2}$
at $R\approx50$ kpc. The median SB probed by KCWI
is about one order of magnitude brighter than that of previous
narrowband survey shown in dashed yellow line in Figure~\ref{fig:SB_r75} and Figure~\ref{fig:SB_r90}.

Let us explore several possibilities that could result in
significantly fainter Ly$\alpha$
profiles that were obtained
by narrowband
surveys at $z\approx2$.
We propose the following explanations:
(1) light loss due to the narrowband filter;
(2) fainter QSOs due to smaller QSO sample constrained by the narrowband.

 The first reason could be due to the Ly$\alpha$ line falling
 outside (or partially outside) the filter.
 In Figure~\ref{fig:2Dspec}, we marked the systemic redshift of these
 QSOs using vertical
 purple lines.
 Similar to \citet{arrigoni19}, for our KCWI sample, we
used the systemic redshift determined from SDSS Mg {\sc ii}\ emission, correcting
the luminosity-dependent offset between Mg {\sc ii}\   and
systemic redshift (Richards et al.
2002; Shen et al. 2016).
Our QSO sample resides at the bright-end in the SDSS QSO sample,
and all of them have a high SNR Mg {\sc ii}\   emission detected in the
SDSS spectra.
 Note that all QSOs except Q1416, have large systemic
 redshift offset comparing to that of the nebular Ly$\alpha$ emission.
More than 50\% of the QSOs have systemic redshifts
$\gtrsim 1000$ km s$^{-1}$  bluer than the Ly$\alpha$
redshift of the extended nebulae.
80\% of the QSOs have $\gtrsim 800$ km s$^{-1}$ velocity offset between
nebular Ly$\alpha$ emisison and QSO systemic redshift determined by SDSS pipeline ($z_{\rm{pipe}}$ in Table~1).
{ If we use the redshift of principle component analysis (PCA) $z_{\rm{PCA}}$ (e.g., \citealt{paris17}), we still
found that 10 out of 16  have velocity offset $>500$ km s$^{-1}$ away from the nebular Ly$\alpha$ redshift. }
When conducting
narrowband surveys, if we require the QSO systemic redshift
to be within 500 km s$^{-1}$ from the most sensitive part of
the filter, then most of nebular emission would be outside the sensitive
part of the filter. For example, the narrowband
filter that Arrigoni Battaia et al. (2016) used has a FWHM of 3000 km s$^{-1}$. Thus,
a velocity offset of $>1300$ km s$^{-1}$ can significantly
  reduce the observed Ly$\alpha$ flux and decrease the ability to
  detect extended emission. { In the Figure~\ref{fig:compare}, we present a simulated
  narrowband image of Q1228}. The systemic redshift of Q1228 is 600 km s$^{-1}$
  bluer than the
  Ly$\alpha$ redshift.
  Using
  a narrowband with FWHM of $\approx 30$ \AA\ centered on the systemic
  redshift (the profile of the filter is the same as the filter used in Arrigoni
Battaia et al. 2016), the Ly$\alpha$ emission is very compact.
  { From Figure~\ref{fig:compare}, if the narrowband filter has a central wavelength
consistent with the systemic redshift (e.g., Arrigoni Battaia et al. 2016),}
Ly$\alpha$ is a compact source which is easily missed under the extremely bright PSF in the narrowband imaging.
The IFU, compared to narrowband,
does not have the narrow redshift constraint, and thus,
one can always detect nebular Ly$\alpha$ emission even if
the  velocity offset between
nebular Ly$\alpha$ and QSO systemic is larger than expected.

Due to the wide wavelength coverage afforded by KCWI,
the KCWI QSO sample is 1.15-mag brighter
than the narrowband QSO sample
(Arrigoni Battaia et al. 2016).
However, using our sample, we demonstrate that there
are no obvious correlation between nebular Ly$\alpha$
and QSO ionizing flux (or QSO magnitude) (see Figure~\ref{fig:LLya_Lnu}). Thus,
there is no strong evidence to support that the 1.15-mag
fainter QSO sample to be a major factor for the one order
of magnitude discrepancy between the extended nebular Ly$\alpha$
observed by KCWI and narrowband (also see the similar
conclusion for $z\sim3$ in \citet{arrigoni19}).

\subsection{Evolution of the Circularly-Averaged Ly$\alpha$ Surface Brightness
from $z\gtrsim3$ to $z\sim2$}

   From Eq.~(1) and Eq.~(3), Figure~\ref{fig:SB_r75} and Figure~\ref{fig:SB_r90},
   we know that 90\% of the $z\approx2$ Ly$\alpha$ nebulae have circularly averaged SB
   profiles fainter than the median SB profiles
   Ly$\alpha$ SB at $z\gtrsim3$ (Borisova et al. 2016; Arrigoni Battaia et al. 2019; also see Marino et al. 2019).
   The median SB at $z\approx2$ is 0.4 dex fainter than the median
   Ly$\alpha$ SB profile at $z\approx3$.
   In this section, we investigate
   possible cause of such an evolution.

 Comparing with $z\approx3$, the lower circularly averaged SB at $z\approx2$
can be arise from two scenarios:
(a) Nebulae at $z\approx2$ could have less circular morphology
or lower emission cover fraction than that at $z\approx3$.
(b) The nebular SB are intrinsically fainter. This could further suggest
 lower local densities or lower gas mass in the QSO halos at $z=2$ compared to $z=3$.
Our measurement is based on a sample of 16 QSOs at $z\approx2$, and
a larger $z\approx2$ QSO sample should be required to further confirm
the results in this paper.

\subsubsection{Comparison between the Covering Fraction of the Ly$\alpha$
Emitting Cloud at $z\gtrsim3$ and $z\approx2$}

For scenario (a), we quantify the
redshift-corrected area of a characteristic
surface brightness. We choose the characteristic
surface brightness of $1\times10^{-17}$ erg s$^{-1}$
cm$^{-2}$ arcsec$^{-2}$ which is well detected in each of
our KCWI QSOs at $z\approx2.3$ and
represented in orange color in Figure~\ref{fig:optimal_extraction}. The fiducial
SB at $z\approx2.3$ is corresponding to $\approx
4\times10^{-18}$ erg s$^{-1}$ cm$^{-2}$ arcsec$^{-2}$
at $z\approx 3.1$. This value is well detected for
each QSOs at $z\approx3.1$ in MUSEUM QSO survey
(Arrigoni Battaia et al. 2019). For all nebulae at $z\approx3$,
the SB contours of
$\gtrsim4\times10^{-18}$ erg s$^{-1}$ cm$^{-2}$ arcsec$^{-2}$
at $z\approx3.1$ have projected size within the KCWI FoV
(Arrigoni Battaia et al. 2019), and thus,
the area at $z\approx2$ and $z\approx3$ can be directly
compared without correcting the difference of the FoV between
KCWI and MUSE.

In Figure~\ref{fig:area}, we present the projected
area of the nebulae with the SB above the fiducial value,
and compared with $z=3$ MUSE results (Arrigoni Battaia et al. 2019)
with the same scaled SB at $z=2$.
From the figure, the median area at $z=2$ is
about 63 arcsec$^2$ (4340 kpc$^2$), this is about 55.1\% of the median
area of MUSE nebulae of $A_{\rm{median}}=123$ arcsec$^2$ (7872 kpc$^2$).
Let us further take into account the halo size from $z\approx3$ -- $z\approx2$.
Let us assume the QSO halo to have a halo mass of $M_{\rm{halo}}\approx10^{12.5}\ M_\odot$
(e.g., White et al. 2012; Shen et al. 2007) at both redshifts.
 Under this assumption, the Virial radius of the halo is $R_{\rm{h, z=3.1}}\approx112$ kpc
 at $z\approx 3.1$, and $R_{\rm{h, z=2.3}}\approx138$ kpc at $z\approx2.3$.
 The {covering} factor of Ly$\alpha$ emitting region
 can be estimated using 4340 kpc$^2$ / ($\pi R_{\rm{h, z=2.3}}$)$\approx0.07$
 at $z\approx2.3$. Using the same method, we estimate that the {covering} factor
 of the Ly$\alpha$ emitting at $z\approx3.1$ is about 0.19, a factor of $2.7\times$ that at $z\approx2.3$.
Therefore, we conclude that QSO halos at $z\approx2$ have a lower
covering factor of the Ly$\alpha$ emitting clouds compared to that at $z\approx3$,
and the lower covering factor may be one of the main reasons for the
evolution of circular-averaged Ly$\alpha$ SB from $z=3$ to $z=2$.

We can further use another independent method to
confirm that the lower circularly-averaged SB at $z\approx2$
seen in Figure~\ref{fig:SB_r75} and Figure~\ref{fig:SB_r90} is due to the lower covering factor.
{We calculated} the median SB in a 90-degree-quadrant for each nebula.
\footnote{Quadrant: a sector of 90 degree, and the quadrant
is chosen with the highest median SB over a large number
of rotations.} Then, we compute the ratio of the 90-degree-quadrant
SB to the full annuli SB. A higher value of this ratio indicates
a lower covering factor or a more asymmetric distribution of the Ly$\alpha$ emitting regions.
In the right panel of Figure~\ref{fig:SB_sector_corr},
we plot the ratio of the
SB in a 90-degree-sector \footnote{Again, for each QSO, the quadrant
is chosen with the highest median SB over a large number
of rotations.} to the median SB in the full annuli.
We could see that this ratio at $z\approx2.3$ (red) is
higher than that at $z\approx3$ (purple) on the scale of 15 -- 70 kpc,
suggesting that the covering factor of the Ly$\alpha$ emitting
clouds at $z\approx2$ is lower than that at $z\approx3$.
In the left panel of Figure~\ref{fig:SB_sector_corr}, we
further show the 90-degree-quadrant, redshift-corrected
SB at $z\approx2.3$. Compared to the SB calculated using full
annuli, the $z\approx2$ and $z\approx3$ SB are more
consistent with each other, indicating that the characteristic, local Ly$\alpha$ SB
within the Ly$\alpha$ emitting regions at $z\approx2$ is close to that at $z\approx3$.
The dimming
of the circularly average SB at $z\approx2$ is due to a smaller
covering factor of Ly$\alpha$ emitting clouds.

\begin{figure}[h]
\includegraphics[width=0.55\textwidth,height=0.32\textheight]{area2.pdf}
 \caption{The distribution of the projected area with redshifted
 surface brightness (SB$_{\rm{z=2.3}}$) greater than $1.0\times10^{-17}$
 erg s$^{-1}$ cm$^{-2}$ arcsec$^{-2}$. Because all MUSE nebulae at $z\approx3$
 have redshifted SB$_{\rm{z=2.3}}\gtrsim 1.0\times^{-17}$ erg s$^{-1}$
 cm$^{-2}$ arcsec$^{-2}$ contour smaller than the KCWI FoV ($16"\times20"$),
 the area at $z\approx2$ and $z\approx3$ can be directly compared without correcting the
 difference of the FoVs between KCWI and MUSE. The median projected area at $z=2.3$
 is 63 arcsec$^2$ which is about 51.3\% of the median
area of $A_{\rm{median}}=123$ arcsec$^2$) for nebulae at $z\approx3.1$. }
 \label{fig:area}
\end{figure}

\begin{figure*}[h!]
\includegraphics[width=1.0\textwidth,height=0.3\textheight]{SB_sector_correction.pdf}
 \caption{{\it Left:} {SB profiles averaged over a 90-degree quadrant chosen as the one with the highest flux}. {\it Right:} the ratio between the
SB in 90-degree-sectors and the median SB in full annuli.
On 15 -- 70 kpc, the ratio at $z\approx2.3$ (red) is
higher than that at $z\approx3$ (purple),
suggesting that the covering factor of the Ly$\alpha$ emitting
clouds at $z\approx2$ is smaller than $z\approx3$ in these
two samples.  }
 \label{fig:SB_sector_corr}
\end{figure*}

\subsubsection{Evolution of the Cosmological Gas Density}

   From the above section, we have shown that at the fixed SB, our KCWI QSO
   sample at $z\approx2$ has less area covered by Ly$\alpha$
 emitting regions compared to QSOs  at $z\approx3$ \citet{borisova16, arrigoni19}.
 In the following, we investigate the scenario (b): whether the $z\approx2$ nebulae
 are intrinsically fainter, i.e., whether $z\approx2$ QSOs have
 lower CGM densities at $z=2$ compared to $z=3$.

 The local surface brightness of the Ly$\alpha$ emitting regions
 between $z\approx2$ and $z\approx3$ is
 consistent with each other (left panel of Figure~\ref{fig:SB_sector_corr}).
 In this section, we further demonstrate that
 cosmic density evolution may not be a major
 factor of causing the evolution of the
 circularly-averaged SB from $z\approx3$ to
 $z\approx2$.

Following the expansion of the Universe,
the cosmic mean density evolves as $(1+z)^{-3}$.
Given the standard convention of defining a dark matter halo
as 200 $\times$ the mean density,
the  density in the halo systematically
decreases from $z=3-2$.
In the standard $\Lambda$CDM Universe, the cosmic scale
factor $a(z)$ is proportional to $\frac{1}{(1+z)}$. At $z>2$, the cosmic
density
$\rho(z)$ is proportional to $\frac{1}{a(z)^3}$, and the density
at $z=3.1$ ($\rho(z=3.1)$) and that at $z=2.3$ ($\rho(z=2.3)$)
have the following relation:
\begin{equation}
\rho(z=3.1)=  \frac{(1+2.3)^3}{(1+3.1)^3}\times \rho(z=2.3) = 1.92\times \rho(z=2.3).
\end{equation}
The average density of a halo at $z=3.1$ is
1.92 $\times$ the halo at $z=2.3$.
If QSO halos follow this general
halo density evolution, and further, if the cool gas
density is proportional to the halo matter density,
then the Ly$\alpha$-emitting cool
gas density in QSO halo at $z=3.1$ is expected
to be 1.92 $\times$ that at $z=2.3$.
Then, the SB$_{\rm{Ly\alpha}}$ could evolve
accordingly as the evolution of cool gas density.

The optically thin photo-ionization model
is a favored scenario
for nebular Ly$\alpha$ emission around AGN (e.g., Heckman
et al. 1991a; Cantalupo et al. 2014;
Arrigoni Battaia et al. 2015, 2016; Cai et al.
2017; see also discussions in Cai et al. 2018).
In the optically thin approximation, the Ly$\alpha$ surface brightness
SB$_{\rm{Ly\alpha}}\propto n_{\rm{H}} N_{\rm{H}}$, where $n_{\rm{H}}$ is the number density of the hydrogen,
and $N_{\rm{H}}$ is the column density which is also
proportional to the gas number density.
From above, we know that the QSO halo density at $z=3.1$
could be 1.92 $\times$ the density of halo at $z=2.3$.
If the gas in the halo is largely optically thin,
then we expect that the $SB_{\rm{Ly\alpha}}\propto n^2$, where $n$ is the gas density. Then, we expect:
\begin{equation}
 SB_{\rm{Ly\alpha}, z=2.3}
= \frac{1}{3.7} \times SB_{\rm{Ly\alpha}, z=3.1},
\end{equation}
The Equ.(8) suggests that if the dimming of the circular-averaged SB is due to the global
evolution of the CGM density, then the characteristic redshift-dimming-corrected
 SB$_{\rm{Ly\alpha}}$ at $z\approx3.1$ could be
a factor of $3.7 \times$ that at $z\approx2.3$.
This is not consistent with our current observations.
As described in the previous section, the
characteristic Ly$\alpha$ emission at $z\approx2$ is  consistent with that at $z\approx3$, but
$z\approx2$ halo has a lower covering factor comparing to $z\approx3$.
This suggests that the local densities of the cool, Ly$\alpha$-emitting gas
have little evolution from $z\approx3$ -- $z\approx2$ as the cosmic expansion,
as expected from cosmological theories (the CGM is not in the linear part of the growth of perturbations).

Note the above conclusion is {drawn} under the optically thin photoionization scenario.
The optically thick photoionization is not favoured by our observations.
If the Ly$\alpha$ nebular emission is mainly contributed by optically thick gas,
then the strong correlation between the luminosity of Ly$\alpha$ ($L_{\rm{Ly\alpha}}$)
and the QSO ionizing luminosities should exist.  Nevertheless, from Figure~\ref{fig:LLya_Lnu}, there are no obvious correlation between
nebular $L_{\rm{Ly\alpha}}$ and the QSO ionizing luminosities in our sample\footnote{The ionizing luminosity of the QSO is
determined by rescaling the composite spectrum
from Lusso et al. (2015) to the
QSO luminosity at the rest-frame 1350 \AA.
If the gas is optically
thick to Lyman continuum photons ($N_{\rm{HI}}$ $\gtrsim  10^{17.2}$ cm$^{-2}$),
the cool gas will behave like a ‘mirror’ and
the SB thus follows the relation $SB_{\rm{Ly\alpha}} \propto
L_{\rm{\nu LL}}$
(Hennawi \& Prochaska 2013).
If the source is bright enough
to keep the gas highly ionized, i.e., the gas is optically thin ($N_{\rm{HI}} <$ 10$^{17.2}$
cm$^{-2}$). It can be shown that (Hennawi \& Prochaska 2013)
the Ly$\alpha$ surface brightness is $SB_{\rm{Ly\alpha}} \propto n_{\rm{H}} N_{\rm{H}}$. Therefore,
in the optically thin regime, the
observed surface brightness should not depend on the luminosity
of the targeted QSOs, but on the density of cool gas.},  indicating optically thick scenario is not favoured by our observations.

 \subsubsection{Evolution of Mass Threshold of Hot Halo}

The decreasing of the Ly$\alpha$ emitting regions from $z\approx3$
to $z\approx2$ may also be interpreted using the mechanisms of the cool
gas penetration in the massive halos \citep{dekel09}.
Dekel et al. (2009) study the penetration of
cold gas into the halo as a function of halo mass
and redshift. They find that, from the simulation, halos at $z\gtrsim3$ with
$M>10^{13}\ M_\odot$ are dominated by hot gas, but at
$z\approx2$, the mass threshold of hot halo decreases to $M\approx10^{12}$ M$_\odot$.
This indicates that QSO halos with a typical halo mass of
$M\sim 10^{12.5}\ M_\odot$
(e.g., White et al. 2012) is more likely to be dominated by hot gas at $z\approx2$ while significant
amount of cool gas
could still penetrate such halos at $z\approx3$. This effect may
reduce the surface brightness of Ly$\alpha$ emission from $z\approx3$ to
$z\approx2$.
Note that the cool gas penetration scenario would
hold for any Ly$\alpha$ powering mechanisms
 (see \citet{arrigoni19}), i.e., the decreasing
of the cool gas penetration in massive halos may yield
a decrease of the Ly$\alpha$ emission.
More quantitative analysis on whether this effect
can decrease the covering factor of cool gas should
be further studied by simulations, combined with
deeper observations on other transition lines.

 \subsection{Ly$\alpha$ morphology and its Evolution
 from $z\gtrsim3$ to $z\approx2$}

\subsubsection{Summary of the Ly$\alpha$ Morphology at $z\approx2$ in our KCWI sample}

 From Borisova et al. (2016) and Arrigoni Battaia et al. (2019), the majority of the
 MUSE radio-quiet nebulae at $z\gtrsim 3$ have symmetric morphology.
 In particular, nebulae with scales smaller than 100 kpc
 have circular morphologies. Nevertheless, our new KCWI observations suggest that nebulae
 at $z\approx2$ have  more irregular
 and asymmetric morphologies. In the following, we briefly comment on the morphology,
 kinematics and surface brightness of  individual nebula.

 From Figure~\ref{fig:optimal_extraction}, we observe that Q2121 has two spatially distinct
 components, with a projected separation of $\approx 85$ kpc.
 From Figure~\ref{fig:2Dspec}, the two components have similar
 velocity and the entire Ly$\alpha$ emission has a FWHM of
 $\approx 500$ km s$^{-1}$. Q0048 and Q1426 both have Ly$\alpha$ emitters
 in the KCWI field of view, at the same redshift with the QSOs,
 suggesting strongly overdense environment in both fields.
 The bright Ly$\alpha$ emitter in  the Q0048 field is 5$''$ west
 of the QSO center. Seven nebulae: Q2121; Q2123;
 Q1227; Q1228; Q1230; Q1416 have projected scales of $>150$ kpc and
 therefore extend beyond the FoV of the KCWI medium slicer. Their spatial extents
 are more similar
 as enormous Ly$\alpha$ nebulae (e.g., Cai et al. 2017; Arrigoni Battaia et al. 2018)
 with the projected size of $\gtrsim200$ kpc.

 Q2127, Q0814, Q0107, Q2125 and Q0848 are
 among the lowest Ly$\alpha$ surface brightness in the sample whose projected
 scales are $<100$ kpc.
 Our KCWI observations
 suggest that nebulae around QSOs at $z\approx2$ may be generally more
 asymmetric, irregular, and having stronger field-to-field variations than that at $z\approx3$.

 \subsubsection{Evolution of the Ly$\alpha$ Morphology}

The morphology of Ly$\alpha$ nebulae contains
information about the cool gas
distribution and geometry in the halos.
 Morphology can also be used to directly compare with cosmological simulations (e.g., Weidlinger et al. 2005;
 Cantalupo et al. 2014). In \S3.4, we quantified the asymmetry of the Ly$\alpha$ light distribution following Arrigoni Battaia et al. (2019).
 From \S3.4, we found that the median
ratio between the
semiminor axis and semimajor axis at $z\approx2$ is
$\alpha_{\rm{z=2}}=0.54$, with a scatter of 0.19. For a comparison, at $z\approx3$,
Arrigoni Battaia et al. (2019) measured the average ratio of  $<\alpha_{\rm{z=3}}>=0.71$, with a scatter of 0.13 (also see Figure~\ref{fig:alpha_size}).
This tentitatively suggests that nebulae at $z\approx2$ may be more asymmetric or clumpy than that at $z\approx3$. Such an asymmetry
could further contribute to the lower average surface brightness at $z\approx2$.
{ This can be further seen from Figure~\ref{fig:2Dspec} which the 2D spectra of the nebulae at
 $z\sim3$ (Borisova) look smooth in the velocity space, while nebulae in the
current survey at $z\sim2$ in our sample are more clumpy.}

Galaxy structure is a powerful method for determining
whether a galaxy is undergoing a recent major merger.
Lotz et al. (2008, 2010) find that the asymmetry
is sensitive to mergers with mass ratios of 1:4 or less.
 The halo morphological merger fraction increases from
 $z\approx3- 2$ and decreases from $z=2-1$.
 The merger fraction
 is highest around $z\approx2$ (e.g., Conselice et al. 2008).
 Lopez-Sanjuan et al. (2009) further point out that the merger
 rate at $z\approx 2$ could be a factor of $2\times$ that at $z\gtrsim3$.
 Such an increased merger fraction at $z\approx2$ may result in
 the asymmetric gas morphology we observe at $z\approx2$.

\subsection{Offsets between the systemic redshift and Ly$\alpha$ redshift}

{ As indicated in \S4.1, we
consider three cases: (1) systemic redshift determined from SDSS pipeline ($z_{\rm{pipe}}$);
(2) systemic redshift determined by Mg {\sc ii}\  emission ($z_{\rm{MgII}}$), correcting
the luminosity-dependent small and known offset between Mg {\sc ii}\ and
systemic redshift (Richards et al.
2002; Shen et al. 2016); and (3) systemic redshift determined by the principal component analysis, with the
reference sample has been chosen to have an automated redshift corresponding to the location of the maximum
of the Mg {\sc ii}\ emission line ($z_{\rm{PCA}}$). These three redshifts are all drawn from the SDSS DR12 QSO catalog
\citep{paris17}. }

The uncertainties of these redshifts
 are known to be $\delta_z \sim 0.003$ ($\approx300$ km s$^{-1}$)
(also see Arrigoni Battaia et al. 2016).
 In Figure~\ref{fig:2Dspec}, we show the systemic redshift determined by SDSS pipeline ($z_{\rm{pipe}}$) (purple
 vertical line) relative to the center of the integrated nebular Ly$\alpha$ emission (marked as zero in the x-axis).
 We further show the center of QSO Ly$\alpha$ emission using the orange vertical line.
 We calculate the velocity shift ($v_{\rm{off}}$) between the Mg {\sc ii}-derived systemic redshift ($z_{\rm{sys}}$) and
 the nebular Ly$\alpha$ redshift ($z_{\rm{Neb, Ly\alpha}}$), and the expression is as follows:
 \begin{equation}
 v_{\rm{off}}= (z_{\rm{sys}}- z_{\rm{Neb, Ly\alpha}})/(1+z_{\rm{Neb, Ly\alpha}})\times c
 \end{equation}
 The velocity shifts range from -3947 to +813 km s$^{-1}$, with a median velocity shift of -1081 km s$^{-1}$.
 In our KCWI sample, only Q1416 and Q2125 have the $z_{\rm{sys}}$
 redshifted than the nebular Ly$\alpha$ redshift ($z_{\rm{Ly\alpha}}$). All other QSOs ($87\%$ of our sample)
 have $z_{\rm{sys}}$ blueshifted than the nebular Ly$\alpha$ redshift.
 Similar results have already been reported in \citet{arrigoni19} at $z\approx3$:
 at $z\approx3$, \citet{arrigoni19} show that the  velocity shift between the
 QSO systemics and the nebular Ly$\alpha$ redshift range
 -6000 km s$^{-1}\le $ $v_{\rm{offset}} \le 2000$ km s$^{-1}$, and
80\% nebulosities have negative velocity
shifts, with the median value of
$\Delta v_{\rm{median}} = -782$ km s$^{-1}$. Similar results and discussions can
also be found in Arrigoni Battaia et al. (2019).
 Note that if we use $z_{\rm{PCA}}$ as the systemic redshift, then the velocity offset ranges from
 -4281 km s$^{-1}$ -- 1463 km s$^{-1}$, with the median velocity offset of -553 km s$^{-1}$, with
 the 1-$\sigma$ scattering of 1417 km s$^{-1}$.
If we use  $z_{\rm{MgII}}$ as the systemic redshift, the velocity shifts range from
-1586 km s$^{-1}$ to $1646$ km s$^{-1}$, with the median offset of -92 km s$^{-1}$,
with the 1-$\sigma$ scattering of 877 km s$^{-1}$.

It is surprising to note that most of the SDSS QSO systemic
redshifts in our sample are
blueshifted compared to
the nebular Ly$\alpha$ emission, {and also, most SDSS QSO systemic
redshift has large offset with respect to the nebular Ly$\alpha$ redshift.
Currently, we still do not
find a good explanation for such a systematic shift.}
We only can conclude that with an intrinsic
uncertainty of 300 km s$^{-1}$
for the empirical calibration of the emission (e.g., Shen et al. 2016), such
systemic redshift estimates may not be optimal for
our ultraluminous QSO sample.
Future efforts are definitely needed to better constrain
this fundamental parameter.
{A better strategy to obtain a more accurate systemic redshift (down to a few tens of km s$^{-1}$)  may require sub-mm observational campaign in future.}
 These sub-mm
observations can help us to evaluate whether the systemic
redshift derived from Mg {\sc ii}\ emission is accurate or not
for such a high luminosity QSO sample at $z\approx2$.

\section{Summary}

In this paper, we conduct a systematic, blind survey of
Ly$\alpha$ nebulae around ultraluminous Type-I QSOs at $z\approx2$
using the Keck Cosmic Web Imager (KCWI).
  {This survey allows us to directly compare with similar integral-field-spectroscopic studies at
$z\gtrsim3$ using MUSE and to study the evolution of the cool gas in massive halos.}
The main conclusion of this paper are as follows:

(1)  We find that {14 out of 16 QSOs}  at $z\approx2$ are
associated with Ly$\alpha$ nebulae with projected linear sizes larger than
50 physical kpc (pkpc). Among them, four nebulae
have large Ly$\alpha$ emission with surface brightness of $SB_{\rm{Ly\alpha}} >10^{-17}$
erg s$^{-1}$ cm$^{-2}$ arcsec$^{-2}$ on $>100$ kpc, and their scales extend beyond the KCWI
field of view.

(2) Our KCWI results suggest that the nebulae at $z\approx2$ are
one order of magnitude brighter and more extended than
previous narrowband surveys indicate (e.g., Arrigoni Battaia et al. 2016).
{ We suggest that this is partially due to the limitation of the narrowband survey technique
and a significant fraction of the diffuse Ly$\alpha$ emission may have be missed in the narrow-band
imaging because of large offsets between the true and the estimated quasar systemic redshift from
the MgII line. (see \S4.1). }
The circularly-averaged Ly$\alpha$ profile is also much more brighter than that
of Lyman break galaxies (LBGs) (e.g., Steidel et al. 2011) and Ly$\alpha$
emitters \cite[e.g.,][]{wisotzki16}. No regular rotational kinematic
patterns have been found in the diffuse Ly$\alpha$ emission for our
QSO sample at $z\approx2$.

(3) We directly measure the circularly-averaged surface
brightness (SB) at $z\approx2$ using the integral field spectroscopy, and we
perform a direct comparison with the $z\approx3$ results using VLT/MUSE \citep{borisova16, arrigoni19}.

 The typical circularly-averaged SB profile of Ly$\alpha$ nebulae around $z\approx2$
 QSOs can be described by a power-law with the slope of -1.8, same as the
 one measured for the nebulae around $z\approx3$ QSOs, however its normalisation
 is about a factor 0.4 dex fainter than the SB profiles at $z\approx3$, after correcting for
 the different redshift dimming.  A larger FoV could further constrain the profile on a larger
 radii, especially to differentiate the exponential and power-law profiles \cite[e.g.,][]{arrigoni19}.

(4) Our analysis suggests that
the Ly$\alpha$ emitting cool gas in the QSO halos at $z\approx2$ may have smaller
cover fraction (see \S4.2). This may be one of the reasons for the lower
circularly average SB from $z\approx3$ to $z\approx2$.
Further, nebulae around $z\approx2$ QSOs appear to have more irregular
and asymmetric morphologies compared to QSO nebulae at
$z\approx3$ as quantified by the ratio between the
semiminor and semimajor axis above the 2$\sigma$
detection level (and denoted by $\alpha$ as in
Arrigoni Battaia et al. 2019).  In particular,
the average $\alpha$ at $z\approx2$ is $\sim40$\%
smaller than that at $z\approx3$.
Taking into account these different morphologies, the
local SB values of the Ly$\alpha$ region, once corrected for redshift-dimming,
become similar between the different redshifts,
especially at the radii $>60$ kpc (see e.g., Figure~\ref{fig:SB_sector_corr}). .

Taken all together, our KCWI results suggest that
from $z\approx 3$ to $z\approx 2$, the covering factor
and possibly the overall mass of cool Ly$\alpha$ emitting CGM
could decrease, assuming the same opening angle at $z\approx2$ and
$z\approx3$. For the Ly$\alpha$ emitting region,
 the typical densities of the cool gas around QSOs do not strongly evolve from
$z\approx3$ to $z\approx2$. Last but not least, we note that there is still
a large scatter and variation in the Ly$\alpha$ surface brightness
and covering factor among individual
systems and that a larger sample would be necessary in order
to obtain a better statistical
analysis and comparison among
different redshifts and QSO properties.

\vspace{0.2in}

{\bf Acknowledgement:}
We acknowledge the discussion of the data reduction and analysis
with Chris Martin, Matt Matuszewski, James D. Neill, Erika Hamden, and Donal Sullivan.
ZC acknowledges the supports provided by NASA through
the Hubble Fellowship grant HST-HF2-51370 awarded by the
Space Telescope Science Institute, which is operated by the
Association of Universities for Research in Astronomy, Inc., for
NASA, under contract NAS 5-26555.
SC gratefully acknowledges support from
Swiss National Science Foundation grant PP00P2\_163824.

\newpage


\part{\partial}
\def\arccosh{\mbox{\rm arccosh}}
\def\tr{\mbox{\rm tr}}
\def\det{\mbox{\rm det}}
\def\sl{\hspace{-5pt}/}
\def\f#1#2{{\textstyle\frac{#1}{#2}}}
\def\ff#1#2#3{{\textstyle\frac{#1}{#2#3}}}
\def\arcsinh{{\rm arcsinh}}
\def\Li{{\rm Li}}
\def\Re{\mbox{\rm Re}}
\def\Ei{{\rm Ei}}
\def\wt{\widetilde}
\def\scs{\scriptstyle}
\thispagestyle{empty}
\newcommand{\mathsym}[1]{{}}
\newcommand{\unicode}[1]{{}}
\newcommand\numberthis{\addtocounter{equation}{1}\tag{\theequation}}
\newcommand{\units}[1]{~\mathrm{#1}}
\begin{document}

\title{Vacuum polarization corrections to the $n=2$ and $n=3$ levels in muonic $^4$He}

\author{Wayne W. Repko}\email{repko@pa.msu.edu}
\affiliation{Department of Physics and Astronomy, Michigan State University, East Lansing, Michigan 48824, USA}
\author{Stanley F. Radford}\email{sradford@brockport.edu}
\affiliation{Department of Physics, The College at Brockport, State
University of New York, Brockport, NY 14420}
\author{Duane A. Dicus}\email{dicus@physics.utexas.edu}
\affiliation{Department of Physics and Center for Particles and Fields,
University of Texas, Austin, TX 78712, USA}

\date{\today}

\begin{abstract}
The $n=2$ and $n=3$ levels for muonic Helium are calculated using a potential that includes all one-loop and recoil effects. Electronic vacuum polarization corrections are calculated using an extension of Kinoshita and Nio method. For $n=2$, the results are $2p_{1/2}-2s_{1/2}=1375.05\pm 1.4$ meV and $2p_{3/2}-2s_{1/2}=1521.65\pm 1.4$ meV, essentially in agreement with the latest summary of the current calculations. The $n=3$ results are summarized in tabular form and give $3p_{1/2}-3s_{1/2}=394.76\pm 0.43$ meV and $3d_{3/2}-3p_{3/2}=111.40$ meV.
\end{abstract}

\maketitle

\section{Introduction}

The tension between the determination of the proton size measured using muonic hydrogen energy levels \cite{nature,science} and electron scattering data \cite{Mohr} has focused renewed attention on accurate calculations of the level splitting in muonic atoms \cite{Borie:1982ax,Borie_1,Pachucki,Pachucki_1,Pohl,Peset}. This discrepancy is also being addressed by using increasingly precise measurements of the electronic energy  levels of ordinary hydrogen \cite{Beyer,Fleurbaey}. Part of the problem in the muonic hydrogen extraction of the proton radius is related to the appropriate modeling of its complicated electromagnetic structure consisting of both electric and magnetic form factors. There are also many more QED corrections since the presence of an additional spin means that there are hyperfine and tensor interactions as well as mixing between energy levels. While these corrections are well understood, they, in some sense, get in the way of figuring out the proper way to model the size parameters of the proton. As a consequence, studying a system without these additional terms could prove helpful in understanding the nuclear corrections.

 The case of muonic $^4$He is considerably simpler  since the nucleus is spinless and there is only one nuclear form factor and a single spin \cite{Borie:2012zz,Diepold}. The underlying potential that describes a muon interacting with a point nucleus is
\begin{eqnarray}\label{V}
V(\vec{r})&=&-Z\al\left\{\frac{1}{r}-\frac{\pi}{2\,m_1^2}\de(\vec{r})- \left(\frac{1}{2\,m_1^2}+\frac{1}{m_1m_2}\right)\frac{\vec{L}\Dot\vec{S}}{r^3} +\frac{1}{2m_1m_2}p_i\left(\frac{\de_{ij}}{r}+\frac{\hat{x}_i\hat{x}_j}{r} \right)p_j\right. \nonumber \\
&&\left.\hspace{30pt}-\frac{\al}{2\,\pi}\left(\frac{1}{m_1^2}+\frac{1}{m_1\,m_2}\right) \frac{\vec{L}\Dot\vec{S}}{r^3}-\frac{4\al}{3}\ln(\mu/\la)\left(\frac{1}{m_1^2}+\frac{2Z}{m_1\,m_2} +\frac{Z^2}{m_2^2}\right)\de(\vec{r})\right. \nonumber \\ [6pt]
&&\left.\hspace{30pt}+\frac{4\al}{3}\left(\frac{1}{m_1^2}\ln(\eta_2)+\frac{Z^2}{m_2^2} \ln(\eta_1)\right) \de(\vec{r})-\frac{2Z\al}{m_1m_2} \left(\frac{m_1^2\ln(\eta_1)-m_2^2\ln(\eta_2)}{m_1^2-m_2^2}\right)\de(\vec{r})\right. \nonumber \\ [6pt]
&&\left.\hspace{30pt}-\frac{7Z\al}{6\pi\,m_1\,m_2}\nabla^2\left(\frac{\ln(\mu\,r+\ga)}{r}\right) \right\}\,.
\end{eqnarray}
The first line in Eq.\,(\ref{V}) is the fine structure contribution to the potential derived from the one photon exchange diagram and is usually called the Breit-Pauli potential. Including the relativistic corrections to the kinetic energy completes the fine structure contribution. The remaining terms are one-loop corrections obtained from the muonium potential of Ref.\cite{GRS} by neglecting the terms containing the muon spin ($\vec{\sigma}_2/2$). Ideally, the spin-independent terms should be verified by calculating the one-loop correction to the Breit-Pauli potential. This potential will be the starting point for the calculation of corrections to both the $n=2$ and $n=3$ levels. The parameters are: $m_1=m_\mu=105.658372$ MeV, $m_2=m_{He}=3728.400128$ MeV and $\al=7.297352566\;10^{-3}$. $\mu$ is the reduced mass, $\la$ is the infra-red cutoff, $\eta_j=m_j/(m_1+m_2)$ and $\ga$ is the Euler constant.

Explicit calculations are given for the $n=2$ case to establish that this approach reproduces the usual results for the $n=2$ Lamb shift. The results for $n=3$ follow in the same way and are summarized in tabular form. In either case, the largest corrections to Eq.\,(\ref{V}) are QED electron vacuum polarization corrections to the Coulomb interaction and modifications due to the nuclear form factor. These will be treated in the next two Sections. The next largest corrections are those from the fine structure portion of the potential and are calculated in Section \ref{FS}. The one-loop corrections are treated in Section \ref{one-loop} and the $n=3$ results are summarized in Section \ref{n=3 lev}. There are also  numerous smaller vacuum polarization corrections to the fine structure terms that are summarized in an Appendix.

\section{Vacuum polarization corrections}

Perhaps the simplest way to treat the vacuum polarization corrections to the Coulomb potential is to use the momentum space approach of Kinoshita and Nio \cite{KN}. They show that the order by order energy corrections can be expressed as
\begin{eqnarray} \label{DelE}
\De E_{n,\ell} &=& -\frac{4\pi\,Z\al}{(2\pi)^3}\int\!\!d^{\,3}k \,\frac{1}{\vec{k}^{\,2}}\left(\Pi^{(2)}_f(\vec{k}^{\,2}) +\Pi^{(2)}_f(\vec{k}^{\,2})\Pi^{(2)}_f(\vec{k}^{\,2})+\Pi^{(4)}_f(\vec{k}^{\,2}) \right.\nonumber \\
&&\left. +\Pi^{(2)}_f(\vec{k}^{\,2})\Pi^{(2)}_f(\vec{k}^{\,2})\Pi^{(2)}_f(\vec{k}^{\,2})+ 2\Pi^{(2)}_f(\vec{k}^{\,2})\Pi^{(4)}_f(\vec{k}^{\,2})+\Pi^{(6)}_f(\vec{k}^{\,2})+ \ldots \right)\ro_{n\ell}(k^2a^2)\,,
\end{eqnarray}
where $a$ denotes the Bohr radius $(\mu\al)^{-1}$ and $\ro_{n\ell}(k^2a^2)$ is the Fourier transform of the non-relativistic probability density
\begin{equation} \label{rho}
\ro_{n\ell}(k^2a^2)=\int\!\!d^{\,3}r\langle |\psi_{n\ell}^m(\vec{r})|^2\rangle e^{-i\vec{k}\cdot\vec{r}}\,.
\end{equation}
The angular brackets indicate an average over the degenerate $m$ values for a given $\ell$. For the $n=2$ states \cite{KN}
\begin{equation}\label{KNrho}
\ro_{2S}(k^2a^2)=\frac{Z^4(Z^4-3k^2a^2 Z^2+2k^4a^4)}
{(Z^2+k^2a^2)^4}\,, \qquad \ro_{2P}(k^2a^2) =\frac{Z^6(Z^2-k^2a^2)}{(Z^2+k^2a^2)^4}\,.
\end{equation}
Note that Eq.\,(\ref{DelE}) contains both reducible e.g., $\Pi^{(2)}_f(\vec{k}^{\,2})\Pi^{(2)}_f(\vec{k}^{\,2})$, and irreducible e.g., $\Pi^{(4)}_f(\vec{k}^{\,2})$, contributions. The sum of these contributions can be obtained if one uses the dispersion representation of K\"allen \cite{GK} and Lehmann \cite{HL}. The electron vacuum polarization corrections to the Coulomb potential take the form
\begin{equation} \label{corr}
D(\vec{k}^{\,2})=-\frac{Ze^2}{\vec{k}^{\,2}}-Ze^2\int_{4m_e^2}^\infty\!\!\frac{d\la}{\la} \frac{\De (\la)}{\la + \vec{k}^{\,2}}\,,
\end{equation}
where
\begin{equation}
\Delta(q^2)=-\frac{(2\pi)^3}{3q^2}\sum_n\de^{(4)}(q-q_n)\langle0|j_\mu(0)|n\rangle \langle n|j^\mu(0)|0\rangle\,.
\end{equation}
For any order in $e^2$, the integral in Eq.\,(\ref{corr}) results in
\begin{equation} \label{Pin}
-Ze^2\int_{4m_e^2}^\infty\!\!\frac{d\la}{\la} \frac{\De^{(n)}(\la)}{\la + \vec{k}^{\,2}}=-\frac{Ze^2}{\vec{k}^{\,2}}\Pi^{(n)}(\vec{k}^{\,2})\,,
\end{equation}
where $\Pi^{(n)}(\vec{k}^{\,2})$ is the sum of all $n/2$ loop vacuum polarization corrections, both reducible and irreducible. In general, the expression for the vacuum polarization correction at order $n$ is
\begin{equation}
\De E_{VP}^{(n)}(n,\ell,Z)=-\frac{4\pi\,Z\al}{(2\pi)^3}\int_{4m_e^2}^\infty \!\!\frac{d\la}{\la} \De^{(n)}(\la)\int\! \frac{d^3k}{\vec{k}^2+\la} \rho_{n\,\ell}(k^2a^2)\,,
\end{equation}
For the $2s$ and $2p$ states, the $d^3k$ integrals are
\begin{equation}
\frac{\pi^2\,Z^3(2a^2\la+Z^2)}{2a(a\sqrt{\la}+Z)^4}\quad {\rm and}\quad \frac{\pi^2\,Z^5}{2a(a\sqrt{\la}+Z)^4}\,.
\end{equation}
With these results, the $2s$ and $2p$ states can be expressed as
\begin{eqnarray}
\De E_{VP}^{(n)}2S &=&  -\frac{\mu\,Z^4\al^2}{4} \int_4^\infty\! \frac{dx}{x}\De^{(n)}(x) \frac{(2\be^2x+Z^2)}{(\be\sqrt{x}+Z)^4}\,, \\
\De E_{VP}^{(n)}2P &=&  -\frac{\mu\,Z^6\al^2}{4} \int_4^\infty\! \frac{dx}{x}\De^{(n)}(x) \frac{1}{(\be\sqrt{x}+Z)^4}\,,
\end{eqnarray}
where $\la=m_e^2\,x$ and $\be=m_e\,a$.
In particular, for the $e\bar{e}$ intermediate state, $\De^{(2)}(x)$ is
\begin{equation}
\De^{(2)}(x)=\frac{\al}{3\pi}(1+2/x)\sqrt{1-4/x}\,\theta(x-4)\,.
\end{equation}
For the fourth order Kallen-Sabry \cite{KS} contribution, $\De^{(2)}(x)$ is replaced by $\De^{(4)}(\sqrt{x}/2)$, where $\De^{(4)}(x)$ is given in Ref.\cite{repko-dicus}. The results are contained in Table \ref{HeTot}. All numerical integrals were evaluated using the NIntegrate routine in Mathematica.

There is also a fourth order contribution to the electron vacuum polarization arising from second order perturbation theory. It can be written as
\begin{equation}
E^{(2)}_{n\ell}=\mu^2Z\al\int_0^\infty\!\!dr\,r^2R_{n\ell}(Zr/a)VP(r)\int_0^\infty\!\! dr'\,r'^2\,g_{n\ell}(Zr/a,Zr'/a)VP(r')R_{n\ell}(Zr'/a)\,,
\end{equation}
where $VP(r)$ is
\begin{equation}
VP(r)=-\frac{Z\al}{r}\int_4^\infty\!\!\frac{dy}{y}\Delta^{(2)}(y)e^{-m_e\sqrt{y}\,r}
\,,
\end{equation}
and $g_{n\ell}(Zr/a,Zr'/a)$ is the radial Coulomb Green's function. For $n=2$ and $\ell=0$ with $x=Zr/a,\,x'=Zr'/a$, $E^{(2)}_{20}$ is
\begin{equation}
E^{(2)}_{20}=\frac{\mu Z^2\al^2}{2}\!\!\int_4^\infty\!\! \frac{dy}{y}\Delta^{(2)}(y)\,\,\int_4^\infty\!\!\frac{dy'}{y'}\Delta^{(2)}(y')\, \,\int_0^\infty\!\!dx\,x(1-x/2)e^{-x/2}e^{-bx}\,\,\int_0^\infty\!\!
dx'\,x'g_{20}(x,x')e^{-b'x'}(1-x'/2)e^{-x'/2}\,.
\end{equation}
Here $b=m_ea\sqrt{y}/Z$ and $b'=m_ea\sqrt{y'}/Z$. The integrals over $x$ and $x'$ can be evaluated analytically and the remaining integrals were evaluated numerically. There is a similar expression for $E^{(2)}_{21}$ with $xe^{-x/2}$ replacing $(1-x/2)e^{-x/2}$ and $24$ replacing $2$. The radial Green's functions $g_{20}(x,x')$ and $g_{21}(x,x')$ are \cite{HP,Hostler,JH,LF,repko-dicus}
\begin{subequations}
\begin{eqnarray}
g_{20}(x,x') &=& Ze^{-(x+x')/2}\left[(2-x)(2-x')\left(\ln(x)+\ln(x')+\frac{(x+x')}{4}+2\ga-\frac{15}{4}-{\rm Ei}(x_<)\right)\right. \nonumber \\
&& \left.+12-2x-2x'-\frac{2}{x}-\frac{2}{x'}+\frac{x}{x'}+\frac{x'}{x}-xx' +(2-x_>)e^{x_<}\left(\frac{1}{x_<}-1\right)\right]\,, \\ [6pt]
g_{21}(x,x') &=& \frac{Z}{3}xx'e^{-(x+x')/2}\left[\ln(x)+\ln(x')+ \frac{(x+x')}{4} +2\ga-\frac{49}{12}-\frac{3}{x}-\frac{3}{x^2}-\frac{2}{x^3} \right.\nonumber \\
&&\left.-\frac{3}{x'}-\frac{3}{x'^{\,2}}-\frac{2}{x'^{\,3}}-{\rm Ei}(x_<) +e^{x_<}\left(\frac{1}{x_<}+\frac{1}{x_<^2} +\frac{2}{x_<^3}\right)\right]\,,
\end{eqnarray}
\end{subequations}
where $x=Zr/a$ and $x'=Zr'/a$. $g_{20}(x,x')$ differs from Ref.\,\cite{LF}. The numerical results are in Table \ref{HeTot}.

Finally, there are relativistic corrections to the one-loop electron vacuum polarization contributions. These were calculated using the results of Ref. \cite{rel}.

\section{Nuclear size corrections}
For a spinless nucleus, the simplest nuclear size connection comes from expanding the charge form factor as
\begin{equation}
F(\vec{k}^2)\sim 1-\frac{1}{6}\langle r^2\rangle\vec{k}^2+\cdots\,,
\end{equation}
where $\langle r^2\rangle\equiv r_N^2$ is the nuclear mean squared charge radius. This leads to an effective potential
\begin{equation}
\De V_N(\vec{r})=\frac{2}{3}\pi Z\al r_N^2\de(\vec{r})\,,
\end{equation}
and an energy correction in the $2s$ state of
\begin{equation}
\De E_N2S=\frac{\mu^3Z^4\al^4 r_N^2}{12}= 105.33\,r_N^2\,{\rm meV\,fm^{-2}}\,,
\end{equation}
for $Z=2$. This result is known to over estimate the correction because any reasonable form factor will contain a term order $r_N^3/a^3$ with the opposite sign. To illustrate this, consider the form factor $F(\vec{k}^2)=(r_N^2\vec{k}^2/12+1)^{-2}$, which is often used to parameterize form factors. The corresponding form of $\De E_N2S$ is
\begin{equation}
\De E_N2S=\frac{\mu^3Z^4\al^4\,r_N^2}{12} \left(1-\frac{5Zr_N}{2\sqrt{3}a}\right)=105.33\,r_N^2(1-0.01097\,r_N)\,.
\end{equation}
The result for helium with $r_N=1.681\pm 0.004$ fm \cite{Diepold} is given in Table \ref{HeTot}. There is very little dependence on the power of the form factor, but the correction is of order $r_N/a$ and its coefficient could be introduced as a parameter. The majority of the error in the nuclear size correction comes from the uncertainty in $r_N$. The value for the nuclear polarizability is taken from Ref. \cite{Sick}.
\begin{center}
\begin{table}[h!]
\begin{tabular}{|l|r|r|r|r|r|}
\hline
  &$2s_{1/2}$ &$2p_{1/2}$ &$2p_{3/2}$ &$2p_{1/2}-2s_{1/2}$ &$2p_{3/2}-2s_{1/2}$  \\
\hline
\hh one loop $m_e$ Vac. Pol.& $-2077.25$ & $-411.46$ & $-411.46$ & $1665.79$
& $1665.79$ \\
\hh one-loop $m_e$ rel. VP  & $-0.85$ & $-0.32$ & $-0.04$ & $0.53$ & $0.81$ \\
\hh Kallen-Sabry Vac. Pol.  & $-15.24$& $-3.67$ & $-3.67$ & $11.57$ & $11.57$ \\
\hh 2nd order $m_e$ Vac. Pol. & $-1.90$ & $-0.19$ & $-0.19$ & $1.71$ & $1.71$ \\
\hh one loop $m_\mu$ Vac. Pol. & $-0.34$ & $0.00$ & $0.00$ & $0.34$ & $0.34$ \\
\hh one loop $m_\pi$ Vac. Pol. & $-0.26$ & $0.00$ & $0.00$ & $0.26$ & $0.26$ \\
\hh &  &  &   &  &  \\
\hh Nuclear size & $292.20\pm 1.4$ & $0.00$ & $0.00$ & $-292.20\pm 1.4$ & $-292.20\pm 1.4$ \\
\hh Nuclear size VP corr. & $2.32$ & $-0.02$ & $-0.02$ & $-2.34$ & $-2.34$ \\
\hh Nuclear polarizability & $-2.35\pm 0.13$ &  &  & $2.35\pm 0.13$ & $2.35\pm 0.13$ \\
\hh &  &  &  &  & \\
\hh Fine structure & $-183.30$ & $-183.00$ & $-37.43$ & $0.29$ & $145.86$ \\
\hh Fine structure VP corr. & $0.86$ & $-0.11$ & $0.02$ & $-0.97$ & $-0.84$ \\
\hh Fine structure VP 2nd order & $0.262$ & $-0.17$ & $-0.02$ & $-0.44$ & $-0.29$ \\
\hh &  &  &  &  &  \\
\hh Order $\al^5$ Spin-Orbit & $0.00$ & $-0.31$ & $0.16$ & $-0.31$ & $0.16$ \\
\hh Order $\al^5$ Lamb shift & $11.57$ & $0.03$ & $0.03$ & $-11.54$ & $-11.54$ \\
\hh & & & & & \\
\hh Total & $-1974.28\pm 1.4$ & $-599.23$ & $-452.64$ & $1375.05\pm 1.4$ & $1521.65\pm 1.4$ \\
\hline
\end{tabular}
\caption{Contributions to the $2p_j-2s_{1/2}$ splitting for muonic helium. The entries are in meV. \label{HeTot}}
\end{table}
\end{center}

\section{Fine structure corrections \label{FS}}

The contributions from the various fine structure terms are
\vspace{6pt}
\begin{eqnarray}
\frac{Z\al\pi}{2\,m_1^2}\langle 2s_{1/2}|\de(\vec{r})|2s_{1/2}\rangle
&=&\frac{\mu^3Z^4\al^4}{16\,m_1^2} \label{fsdelta} \\ [6pt]
Z\al\left(\frac{1}{2\,m_1^2}+\frac{1}{m_1m_2}\right)\langle 2p_j| \frac{\vec{L}\Dot\vec{S}}{r^3}|2p_j\rangle &=& \frac{\mu^3Z^4\al^4}{48} \left(\frac{1}{2\,m_1^2}+\frac{1}{m_1m_2}\right)\left(j(j+1)-11/4\right)\label{fsls} \\ [6pt]
-\frac{Z\al}{2m_1m_2}\langle 2s|p_i\left(\frac{\de_{ij}}{r} +\frac{\hat{x}_i\hat{x}_j}{r} \right)p_j |2s\rangle &=& -\frac{3\mu^3Z^4\al^4}{16m_1m_2} \\ [6pt]
-\frac{Z\al}{2m_1m_2}\langle 2p|p_i\left(\frac{\de_{ij}}{r} +\frac{\hat{x}_i\hat{x}_j}{r} \right)p_j |2p\rangle &=& -\frac{\mu^3Z^4\al^4}{16m_1m_2} \\ [6pt]
-\frac{1}{8}\left(\frac{1}{m_1^3}+\frac{1}{m_2^3}\right)\langle 2s| (\vec{p}^{\,2})^2|2s\rangle &=& -\frac{13}{128}\mu^4Z^4\al^4 \left(\frac{1}{m_1^3}+\frac{1}{m_2^3}\right) \\ [6pt]
-\frac{1}{8}\left(\frac{1}{m_1^3}+\frac{1}{m_2^3}\right)\langle 2p| (\vec{p}^{\,2})^2|2p\rangle &=& -\frac{7}{384}\mu^4Z^4\al^4 \left(\frac{1}{m_1^3}+\frac{1}{m_2^3}\right)
\end{eqnarray}
\newpage
Using these results, the total fine structure contributions are
\begin{eqnarray}
E_{FS}(2s_{1/2}) &=& -\frac{13}{128}\mu^4Z^4\al^4 \left(\frac{1}{m_1^3}+\frac{1}{m_2^3}\right)-\frac{3\mu^3Z^4\al^4}{16m_1m_2} +\frac{\mu^3Z^4\al^4}{16\,m_1^2} \\ [4pt]
E_{FS}(2p_{1/2}) &=&  -\frac{7}{384}\mu^4Z^4\al^4 \left(\frac{1}{m_1^3}+\frac{1}{m_2^3}\right)-\frac{\mu^3Z^4\al^4}{24} \left(\frac{1}{2\,m_1^2}+\frac{1}{m_1m_2}\right)-\frac{\mu^3Z^4\al^4}{16m_1m_2} \\ [4pt]
E_{FS}(2p_{3/2}) &=& -\frac{7}{384}\mu^4Z^4\al^4 \left(\frac{1}{m_1^3}+\frac{1}{m_2^3}\right)+\frac{\mu^3Z^4\al^4}{48} \left(\frac{1}{2\,m_1^2}+\frac{1}{m_1m_2}\right)-\frac{\mu^3Z^4\al^4}{16m_1m_2}
\end{eqnarray}
The numerical results are given in Table \ref{HeTot}. Note that the $2p_{1/2}$ and $2s_{1/2}$ are not degenerate due to recoil effects.

\section{One-loop contributions \label{one-loop}}
The one-loop corrections consist of a spin-orbit contribution and a spin-independent contribution, both of order $Z^4\al^5$. The spin-orbit correction is
\begin{equation}
\langle 2p_j|Z\al a_\mu\left(\frac{1}{m_1^2}+\frac{1}{m_1m_2}\right) \frac{\vec{L}\Dot\vec{S}}{r^3}|2p_j\rangle = \frac{\mu^3Z^4\al^4a_\mu}{48} \left(\frac{1}{m_1^2}+\frac{1}{m_1m_2}\right)(j(j+1)-11/4)\,,
\end{equation}
where $a_\mu$ is the muon anomalous magnetic moment.

The matrix element of the remaining one-loop spin independent term $V'_{SI}$ is \cite{GRS}
\begin{eqnarray}
\langle n\ell|V'_{SI}|n\ell\rangle &=& \frac{\mu^3Z^4\al^5}{n^3\pi} \left\{\frac{4}{3}\left(\left(\log\left[\frac{1}{Z^2\al^2k_0[n,0]}\right]+\frac{5}{6}\right) \left(\frac{1}{m_1^2}+\frac{2Z}{m_1m_2}+\frac{Z^2}{m_2^2}\right)\right.\right. \nonumber \\
&&\left.\left. -\left(\frac{1}{m_1^2}\log(\eta_2)+\frac{Z^2}{m_2^2}\log(\eta_1)\right)\right)+\frac{14Z}{3m_1m_2} \Big(\log\left[\frac{2Z\al}{n}\right]+\frac{n-1}{2n}+\sum_{j=0}^n\frac{1}{j}\Big)\right. \nonumber \\
&&\left.+\frac{2Z}{m_1m_2}\left(\frac{m_1^2\log(\eta_1)-m_2^2\log(\eta_2)}{(m_1^2-m_2^2)}\right) \right\}\de_{\ell 0} \nonumber \\
&&+\frac{\mu^3Z^4\al^5}{n^3\pi} \left\{\frac{4}{3}\log\left[\frac{1}{k_0[n,\ell]}\right]\left(\frac{1}{m_1^2}+\frac{2Z}{m_1m_2} +\frac{Z^2}{m_2^2}\right)-\frac{7Z}{3m_1m_2}\frac{1}{\ell(\ell+1)(2\ell+1)}\right\} (1-\de_{\ell 0})
\end{eqnarray}

The Bethe sums $k_0[n,\ell]$ used are $\ln(k_0[1,0])=2.9841285$, $\ln(k_0[2,0])=2.8117699$ and $\ln(k_0[2,1])=-0.0300167$ \cite{KM}.
Again, the numerical results are listed in Table \ref{HeTot}.

\section{Results for $n=3$ \label{n=3 lev}}

The corresponding results for $n=3$ can be calculated using the same techniques that were used in the $n=2$ case. Wave functions for the relativistic corrections can be derived from Rose \cite{Rose}, Eqs.\,(5.49a) and (5.49b). The necessary radial Green's functions are \cite{LF}
\begin{subequations}
\begin{eqnarray}
g_{30}(x,x')&=&\frac{4Z}{9}e^{-(x+x')/2}\left[\left(x^2/2-3x+3\right)\left(x'^2/2-3x'+3\right) \left(\log(x)+\log(x')+(x+x')/6+2\ga-25/6-\Ei(x_<)\right)\right.\nonumber \\
&& \left.+\left((3-x)x/3-5x/2-1/x+7\right)\left(x'^2/2-3x'+3\right)+\left(x^2/2-3x+3\right) \left((3-x')x'/3-5x'/2-1/x'+7\right)\right. \nonumber \\ [4pt]
&& \left. +\left(x_>^2/2-3x_>+3\right) e^{x_<}\left(x_</2+1/x_<-5/2\right)\right]\,, \\
g_{31}(x,x')&=&\frac{Z}{18}xx'e^{-(x+x')/2}\left[(4-x)(4-x') \left(\log(x)+\log(x')+(x+x')/6+2\ga-55/12-\Ei(x_<)\right)\right.\nonumber \\
&&\left.(4-x)\left(x'/3-6/x'-4/x'^2-2/x'^3+5\right)+ \left(x/3-6/x-4/x^2-2/x^3+5\right)(4-x')\right. \nonumber \\ [4pt]
&&\left.+(4-x_>)e^{x_<}\left(2/x_<^3+2/x_<^2+3/x_<-1\right)\right]\,, \\
g_{32}(x,x')&=&\frac{Z}{90}x^2x'^2e^{-(x+x')/2}\left[\log(x)+\log(x')+(x+x')/6+2\ga-89/20-\Ei(x_<) \right.\nonumber \\
&& \left.-5/x-10/x^2-20/x^3-30/x^4-24/x^5-5/x'-10/x'^2-20/x'^3-30/x'^4-24/x'^5 \right.\nonumber \\
&&\left. +e^{x_<}\left(1/x_<+1/x_<^2+2/x_<^3+6/x_<^4+24/x_<^5\right)\right]\,,
\end{eqnarray}
\end{subequations}
where $x=2Zr/3a$ and $x'=2Zr'/3a$. The $n=3$ Bethe sums can be found in Ref.\cite{KM}

The resulting level corrections and splittings are summarized in Table \ref{n=3}.
\begin{table}[h]
\centering
\begin{tabular}{|l|r|r|r|r|r|r|r|}
\hline
  &$3s_{1/2}$ &$3p_{1/2}$ &$3p_{3/2}$ &$3d_{3/2}$ &$3d_{5/2}$ &$3p_{1/2}-3s_{1/2}$ &$3d_{3/2}-3p_{3/2}$  \\
\hline
\hh one loop $m_e$ Vac. Pol.& $-603.07$ & $-121.83$ & $-121.83$ &$-11.37$ &$-11.37$ &$481.24$ &$110.46$ \\
\hh one-loop $m_e$ rel. VP  & $-0.26$ &$-0.10$ &$-0.02$ &$0.00$ &$0.00$ & $0.16$ &$0.02$ \\
\hh Kallen-Sabry Vac. Pol.  & $-4.41$ &$-1.05$ & $-1.05$ &$-0.13$ &$-0.13$ & $3.36$ &$0.92$ \\
\hh 2nd order $m_e$ Vac. Pol. & $-0.86$ &$-0.10$ &$-0.10$ &$0.00$ &$0.00$ &$0.76$ &$0.10$ \\
\hh one loop $m_\mu$ Vac. Pol. &$-0.10$ &$0.00$ &$0.00$ &$0.00$ &$0.00$ & $0.10$ & $0.00$ \\
\hh one loop $m_\pi$ Vac. Pol. &$-0.08$ &$0.00$ &$0.00$ &$0.00$ &$0.00$ & $0.08$ & $0.00$ \\
\hh &  &  &   &  &  &  &  \\
\hh Nuclear size &$86.58\pm 0.41$ &$0.00$ &$0.00$ &$0.00$ &$0.00$ & $-86.58\pm 0.41$ & $0.00$ \\
\hh Nuclear size VP corr. & $1.04$ &$-0.01$ &$-0.01$ &$0.00$ &$0.00$ & $-1.05$ & $0.01$ \\
\hh Nuclear polarizability &$-0.70\pm 0.13$  & 0.00 & 0.00 & 0.00 &0.00  & $0.70\pm 0.13$ & $ $ \\
\hh &  &  &  &  &  &  &  \\
\hh Fine structure & $-65.01$ &$-64.92$ & $-21.79$ &$-21.84$ &$-7.39$ & $0.09$ & $-0.05$ \\
\hh Fine structure VP corr. & $0.25$ &$-0.06$ &$-0.01$ &$0.00$ &$0.00$ & $-0.31$ &$0.01$ \\
\hh Fine structure VP 2nd order & $0.20$ &$-0.09$ &$-0.01$ &$0.00$ &$0.00$ & $-0.29$ & $0.01$ \\
\hh &  &  &  &  &  &  & \\
\hh Order $\al^5$ Spin-Orbit & $0.00$ & $-0.07$ &$0.03$ &$-0.02$ &$0.01$ & $-0.07$ & $-0.05$ \\
\hh Order $\al^5$ Lamb shift & $3.45$ &$0.01$ &$0.01$ &$0.00$ &$0.00$ & $-3.44$ &$-0.01$ \\
\hh & & & & & &  &\\
\hh Total & $-582.97\pm 0.43$ &$-188.20$ &$-144.76$ &$-33.37$ &$-18.88$ &$394.76\pm 0.43$ & $111.40$ \\
\hline
\end{tabular}
\caption{Contributions to the $3p_{1/2}-3s_{1/2}$ and $3d_{3/2}-3p_{3/2}$ splittings for muonic helium. The entries are in meV. \label{n=3}}
\end{table}
The nuclear polarizability correction was obtained by rescaling the result in Ref.\cite{Sick} using the $3s$ wave function evaluated at the origin.

\section{Conclusions}

Starting from a potential describing the QED interaction between a point scalar nucleus and a muon that includes the one-loop level corrections and all recoil effects, the Lamb shifts for $n=2$ and $n=3$ muonic helium have been calculated. The effects of electron vacuum polarization were included by extending the Kinoshita-Nio \cite{KN} method to all operators in the potential. This significantly simplifies the calculation of the various corrections. The nuclear size correction was modeled using a simple one-parameter form factor. Given an accurate value of $r_N$, the resulting correction can be included in the prediction of the $2p_{1/2}-2s_{1/2}$ and $2p_{3/2}-2s_{1/2}$ splittings and the corresponding results for $n=3$, or the correction can be used to extract a value of $r_N$ from spectroscopic data.