# Extragalactic CO emission lines in the CMB experiments: a forgotten signal
and a foreground

Abhishek S. Maniyar<EMAIL_ADDRESS>Kavli Institute for Particle
Astrophysics and Cosmology, Stanford University, 452 Lomita Mall, Stanford, CA
94305, USA  SLAC National Accelerator Laboratory, 2575 Sand Hill Road, Menlo
Park, CA 94025, USA  Center for Cosmology and Particle Physics, Department of
Physics, New York University, New York, NY 10003, USA  Athanasia Gkogkou  Aix
Marseille Univ, CNRS, CNES, LAM, Marseille, France  William R. Coulton  Center
for Computational Astrophysics, Flatiron Institute, 162 5th Avenue, New York,
NY 10010, USA  Zack Li  Canadian Institute for Theoretical Astrophysics,
University of Toronto, Toronto, ON, Canada M5S 3H8  Guilaine Lagache  Aix
Marseille Univ, CNRS, CNES, LAM, Marseille, France  Anthony R. Pullen  Center
for Cosmology and Particle Physics, Department of Physics, New York
University, New York, NY 10003, USA

###### Abstract

High resolution cosmic microwave background (CMB) experiments have allowed us
to precisely measure the CMB temperature power spectrum down to very small
scales (multipole $\ell\sim 3000$). Such measurements at multiple frequencies
enable separating the primary CMB anisotropies with other signals like CMB
lensing, thermal and kinematic Sunyaev-Zel’dovich effects (tSZ and kSZ), and
cosmic infrared background (CIB). In this paper, we explore another signal of
interest at these frequencies that should be present in the CMB maps:
extragalactic CO molecular rotational line emissions, which are the most
widely used tracers of molecular gas in the line intensity mapping
experiments. Using the SIDES simulations adopted for top hat bandpasses at 150
and 220 GHz, we show that the cross-correlation of the CIB with CO lines has a
contribution similar to the CIB-tSZ correlation and the kSZ power, thereby
contributing a non-negligible amount to the total power at these scales. This
signal, therefore, may significantly impact the recently reported $\geq
3\sigma$ detection of the kSZ power spectrum from the South Pole Telescope
(SPT) collaboration, as the contribution of the CO lines is not considered in
such analyses. Our results also provide a new way of measuring the CO power
spectrum in cross-correlation with the CIB. Finally, these results show that
the CO emissions present in the CMB maps will have to be accounted for in all
the CMB auto-power spectrum and cross-correlation studies involving a LSS
tracer.

##  I Introduction

Observations of the cosmic microwave background (CMB) continue to help us
understand the properties of the Universe. A distinct advantage of CMB
observations is that different CMB anisotropies contain information coming
from different epochs of the Universe. For example, the primary CMB
anisotropies which show up on large angular scales in the CMB power spectrum
($\ell<2000$) precisely constrain properties of the Universe around redshift
$z\sim 1100$ Planck_cosmo_2020 . On the other hand, secondary anisotropies
like CMB lensing, thermal and kinematic Sunyaev-Zel’dovich effects (tSZ and
kSZ) are tracers of the low redshift Universe. They show up on relatively
smaller angular scales ($\ell>2000$) and, as a result, are hard to measure and
disentangle. Additional sources of emission, such as the cosmic infrared
background (CIB), sometimes dominate over these anisotropies and make these
measurements even more difficult. However, due to the wealth of crucial
information they contain, precisely measuring these small-scale CMB
anisotropies is a major goal of current and future CMB experiments. This
requires careful modeling of all sources of anisotropies and emissions present
at such scales.

Fortunately, the ever-increasing sensitivity and angular resolution of CMB
experiments like Planck Plancktsz_2016 , ACT Dunkley_2013 , SPT George_2015 ,
SO Ade_19 , CMB-S4 Abazajian_16  have brought us closer to this goal. Although
foregrounds like the CIB contaminate CMB maps near the frequencies where the
CMB dominates ($\sim 150$ GHz), they have a distinct spectral shape with
respect to the CMB. As a result, an effective way to constrain different
components of the CMB map is to perform the measurement across different
frequencies and utilize this distinct spectral dependence. These experiments
have done exactly this by using their multi-frequency data to successfully
separate different individual contributions to the CMB power spectrum at these
small scales. In fact, for the first time, the latest results from the SPT
experiment measured the evasive kSZ power spectrum at $\geq 3\sigma$
Reichardt_2021 . Measuring the kSZ power spectrum is challenging due to its
small amplitude and the presence of other sources like the tSZ, and CIB
Dunkley_2013 ; George_2015 ; Reichardt_2021 .

The kSZ signal contains crucial information about the reionization history of
the Universe, which is otherwise extremely difficult to obtain. The SPT’s
achievement, therefore, has paved a way for future, more sensitive
measurements of the kSZ power spectrum and, as a result, a better
understanding of the epoch of reionization. In light of the upcoming more
sensitive CMB data, it is imperative to carefully scrutinize the current
analysis methods to perform a comprehensive study by considering all possible
sources of power at these frequencies. Reichardt_2021  perform a thorough
analysis of the small-scale SPT data at 95, 150, and 220 GHz. Apart from the
primary CMB, they consider the following templates and fit for them: tSZ, CIB,
CIB-tSZ cross-correlation, radio-source emission, galactic dust, and kSZ.

In this paper, we explore another signal of interest at these frequencies:
bright CO molecular rotational line emissions. They have garnered interest
especially in the field of line intensity mapping (LIM). LIM experiments aim
to make a statistical observation of the aggregate emission from many
unresolved galaxies rather than aiming to image individual galaxies directly
(see Kovetz_2017 ; Bernal_2022  for a review). In LIM studies, CO lines are
the most widely used tracers of molecular gas, which is a crucial interstellar
medium component. These rotational emission lines are among the brightest in
the galaxy spectra and occur at the ladder of frequencies $\nu_{J\rightarrow
J-1}=J\times 115.27$ GHz for the $J\rightarrow J-1$ transitions. Different CO
transitions from different redshifts will be observed in CMB bands ($\sim 150$
GHz). As first pointed out by Righi_2008 , the CO lines should be a source of
foreground for the CMB observations. However, they also point out that the
amplitude of the CO power spectrum goes down with the spectral resolution of
the instrument i.e. the amplitude goes down with lower values of
$\Delta\nu/\nu_{\rm obs}$, where $\Delta\nu$ is the bandwidth of the
instrument at the observed frequency $\nu_{\rm obs}$. As a result, for a SPT-
like instrument with $\Delta\nu\approx 50$ GHz, the amplitude of different CO
lines acting as a foreground should be small Righi_2008 . For example,
according to the model used by Righi_2008  for different CO lines, with
$\Delta\nu/\nu_{\rm obs}=0.2$ at $70$ GHz, all the CO lines have their
respective power spectrum amplitudes
$\mathcal{D}_{\ell}=\ell(\ell+1)C_{\ell}/2\pi\><0.1\,\mu{\rm K}^{2}$ at
$\ell=3000$. By comparison, the smallest amplitude signal of interest in the
CMB power spectrum is the kSZ signal, and it is expected to be of the order
$\mathcal{D}_{\ell}\sim 1\mu{\rm K}^{2}$ around $\ell=3000$ George_2015 ;
Reichardt_2021 . Perhaps, for this reason, the contribution of the CO lines to
the CMB power spectrum has been neglected so far in the analysis.

However, as shown in Padmanabhan_2018 ; Keating_2016 , predictions for CO line
power spectra from different models differ by more than an order of magnitude.
This is mainly due to our limited knowledge of the astrophysics of the
galaxies. The main goal of this paper is to revisit this signal and compare
the amount of CO line emissions in the CMB maps to other signals/foregrounds
at small scales. Using the “simulated infrared dusty extragalactic sky”
(SIDES) simulation Bethermin_2017 ; Bethermin_2022 , we show that CO lines
have a non-zero correlation with other signals like the CIB and tSZ. Our most
important point is that, already at the level of the current CMB data, the
total CO signal in combination with these signals is significant. In fact, we
show that the sum of the cross-power spectra between the CIB and different CO
lines can have an amplitude similar to that of the currently estimated kSZ
power spectrum level. This result, therefore, paves a new way to measure the
CO signal from the CMB data. It may also have a major impact on the potential
measurement of the kSZ power spectrum. Our result thus motivates including
such signals in all future small-scale CMB data analyses.

The paper is organized as follows. We first present in Sect. II the SIDES
simulations used to produce, in a consistent way and for the first time, the
CO (total and individual lines), CIB, and tSZ maps. The maps are then used to
compute the auto- and cross-spectra of these components, which are discussed
in Sect. III. This section reveals that the CIB-CO correlation is a bright
foreground component i) that has to be considered to revisit the kSZ
measurements (Sect. IV.1), ii) that can be used to probe the CO signal (Sect.
IV.2) and iii) that has to be included in any cross-correlation analysis
(Sect. IV.3). The conclusion is given in Sect. V. We also discuss in Appendix
A the contribution of the different rotational CO lines to their cross-
correlation with tSZ and CIB.

##  II Maps from SIDES Simulations

SIDES (Simulated Infrared Dusty Extragalactic Sky) is a publicly
available111https://gitlab.lam.fr/mbethermin/sides-public-release simulation
of the far-infrared (FIR) and sub-millimeter sky based on observed empirical
relations Bethermin_2017 ; Bethermin_2022 . A detailed description of the
simulation is given in Bethermin_2017 ; Bethermin_2022 ; gkogkou2022 . We
provide a brief description of the model in the following sections.

###  II.1 Continuum emission from dusty star-forming galaxies in SIDES

The starting point of SIDES is a dark matter lightcone, which is given by the
Uchuu N-body simulation ishiyama2021 . This cosmological simulation achieves
both a high mass resolution ($\rm 3.27\times 10^{8}\,M_{\odot}h^{-1}$) and a
large comoving volume (with a comoving side-length box of 2000 $\rm
Mpc\>h^{-1}$) between $0<z<7$. The dark matter halos are populated with
galaxies of certain stellar mass through abundance matching. The generated
galaxies are then split into passive and star-forming with a probability
determined by observations davidzon2018 . It is assumed that only the star-
forming galaxies emit in the FIR and millimeter, so only these types of
galaxies need to be assigned with a star formation rate (SFR) value. The SFR
values of the main sequence and starburst galaxies are drawn accordingly based
on the parameterized fit of the SFR-$M_{\rm stellar}$ relation described in
schreiber2015  taking into account a scatter of 0.3 dex. The $L_{\rm IR}$ of
each galaxy is subsequently defined by the drawn SFR value using the Kennicutt
conversion factor ($1.0\times 10^{-10}M_{\odot}L_{\odot}^{-1}{\rm yr}^{-1}$)
kennicutt1998 . The $L_{\rm IR}$ normalizes the spectral energy distribution
(SED) of the galaxy. The shape of the SED depends on the galaxy type (main
sequence or starburst) and on the mean intensity of the radiation field
($\langle U\rangle$), which is correlated to the temperature of the dust
(magdis2012, ; bethermin2015, ; Bethermin_2017, ).

###  II.2 CO lines in SIDES

SIDES also simulates the emission of the strongest high-redshift
(sub-)millimeter lines such as [CII], CO, and [CI] (Bethermin_2022, ). In this
work, we focus only on the CO molecular lines. However, other high redshift
lines like [CI] are potentially present in the CMB maps as well and should be
studied similarly. For CO, the fundamental transition is modeled in SIDES from
the $L_{\rm IR}-L_{\rm CO(1-0)}^{{}^{\prime}}$ correlation (sargent2014, ) for
the main sequence galaxies, while for the starburst systems there is an offset
of -0.46 dex for $L_{\rm IR}$ at a given $L_{\rm CO(1-0)}^{{}^{\prime}}$. The
flux of the other CO transitions is computed using a clumpy and diffuse
spectral line energy distribution (SLED) template from bournaud2015  and
following the empirical relation presented in daddi2015  which connects the
flux ratio of CO(5-4) and CO(2-1) transitions with the mean intensity of the
radiation field $\langle U\rangle$.

###  II.3 Astrophysical models for tSZ from SIDES halos

We apply a semi-analytic prescription for the thermal Sunyaev-Zel’dovich (tSZ)
effect closely following websky , in which we paste gas pressure profiles
obtained from detailed hydrodynamic simulations onto the halo catalog obtained
from the Uchuu N-body simulation (i.e. as in sehgal2009, ; trac2011, ;
osatonagai2022, ). The tSZ induces a spectral distortion in CMB photons
proportional to the line-of-sight integral of the ratio of electron thermal
and rest energies,

$\Delta T(\nu)=T_{\mathrm{CMB}}\,f(x)\int
n_{e}\frac{k_{B}T_{e}}{m_{e}c^{2}}\sigma_{T}dl\equiv
T_{\mathrm{CMB}}\,f(x)\,y.$ (1)

Here $x=h\nu/k_{B}T_{\mathrm{CMB}}$ and the frequency dependence relative to
the CMB blackbody in the non-relativistic limit is

$f(x)=\left(x\frac{e^{x}+1}{e^{x}-1}-4\right).$ (2)

The evolution of cluster gas is governed by the physics of star formation and
AGN feedback, and as in websky  we include these effects through a parametric
model for the radial cluster pressure profile based on halo mass and redshift,
calibrated from detailed hydrodynamic simulations (BBPS2012, ). We apply the
profile with parameters from Table 1 of BBPS2012  corresponding to AGN
feedback $\Delta=200$. These are implemented in a publicly available script
using the package XGPaint.jl222https://github.com/xzackli/tsz_sides.

###  II.4 Why use SIDES?

The SIDES model is in very good agreement with all the observations we can
find for dusty star-forming galaxies (DSFG) in the whole redshift range of
SIDES ($0<z<7$). In addition, the model also reproduces more general
quantities, such as the stellar mass function (SMF), the evolution of sSFR
Bethermin_2017 . For the purpose of this paper (CO and CIB), the specific
observations SIDES reproduces include:

  * •

The CO line luminosity functions that are observed in the redshift range
$1<z<3.4$ (Fig. 5 and 6 of gkogkou2022 ).

  * •

The measurement of the shot noise level of the CO emission at 3 mm (Fig. 18 of
Bethermin_2022 ).

  * •

The power spectra of the CIB anisotropies (Fig. 3 and 4 of gkogkou2022 ).

  * •

The number counts from Spitzer, Herschel and (sub-)millimeter experiments
(e.g. SCUBA2, ALMA), galaxy redshift distributions, number counts per redshift
slice and P(D) from SPIRE/Herschel Bethermin_2017 .

  * •

The obscured star formation rate density (SFRD) measurements up to z=4 (Fig. 3
of Bethermin_2022 ). Notice that the observational constraints at $z>4$ are
very dispersed.

For tSZ, as discussed in the next section, the power spectrum from our
simulation is compatible (within 1$\sigma$) with the recent measurements.

Thus SIDES produces consistent CIB, CO maps, and tSZ maps, both for their
emission signals and clustering properties on $\sim$120 square degrees. It is
then easy to compute the cross-power spectra to estimate their contribution to
the small-scale CMB anisotropy measurements. These are the first such
simulations consistently containing all these components.

##  III Results

###  III.1  Comparison between different foreground power spectra

In this work, we simulate all the components as they would be observed using
the SPT telescope. For the CIB emission, we bandpass the galaxy SEDs with SPT-
SZ filters at 150 and 220 GHz. On the other hand, to go from frequency
independent y map to the tSZ emission corresponding to the 150 and 220 GHz
SPT-SZ bandpass, we multiply it by the numbers provided in Tab. A3 of
Lagache_2020  (which are -0.416 and 9.44 respectively). To generate the CO
emission maps, we use an approximate (but fast) method. We assume that the
SPT-SZ bandpass at 150 and 220 GHz has uniform spectral response between
125-175 and 195-245 GHz respectively. With this, we then simulate the CO map
centered at 150 and 220 GHz with a spectral bandwidth of $\Delta\nu=50$ GHz.
We checked on a 2 Sq. Deg. simulated patch using the measured 150 GHz SPT-SZ
bandpass that this approximation has no impact on the CO power spectra. In
this section, we present the power spectra of these signal maps and their
cross-correlations. It is important to note in our model that the CIB maps
include all components of dusty star-forming galaxies, i.e. 1- halo and 2-halo
terms (usually called CIB), and shot noise. We also do not implement any flux
cut in our CIB maps.

Figure 1:  tSZ (y) power spectrum from SIDES catalog (black) along with data
from Planck Plancktsz_2016 ; Bolliet_2018  (blue and orange), SPT
Reichardt_2021  (green), and ACT Choi_2020  (red) experiments after removing
40 halos massive than $5\times 10^{14}M_{\odot}$. With this choice, the power
spectrum from our simulation is compatible with these measurements for
$\ell>400$, i.e. in the range of scales of interest. If these halos are not
removed, tSZ power spectrum amplitude from our simulations is much higher than
the observations. This is due to the sample variance caused by the small
volume of our simulations.  Figure 2: Power spectrum (absolute value) of
different components of the SIDES simulation for SPT 150 and 220 GHz channels.
The curves for kSZ and radio components have been taken from Reichardt_2021 .
The CIB-labeled curves include CIB and Poisson. As can be seen, CIB-CO cross-
correlation is quite significant at both frequencies. In fact, at 220 GHz,
CIB-CO correlation is the second most important signal after the CIB.

First, we verify that our tSZ profile painting method on the dark matter halos
works correctly. For this, we plot the frequency-independent power spectrum of
our tSZ maps and compare it with the corresponding Planck Plancktsz_2016 ;
Bolliet_2018 , SPT Reichardt_2021 , and ACT Choi_2020  data in Fig. 1 in
dashed black curve. We can see that on the relevant small scales here, the tSZ
power spectrum from our maps is higher than the measured values. Due to its
mass weighting, the tSZ power spectrum is very sensitive to Poisson
fluctuations Komatsu_2002 . The high signal seen here is because of the sample
variance due to the small volume of our simulation ($\sim 120$ square
degrees). After removing halos more massive than $5\times 10^{14}M_{\odot}$,
we get a tSZ power spectrum more compatible with the observations (solid black
line) 333The code used to produce this tSZ map has been extensively tested
including reproducing the Websky Stein_2020  simulations.. Therefore, we use
this catalog with 40 halos $5\times 10^{14}M_{\odot}$ removed. We have checked
that this choice has insignificant impact on the CIB and CO power spectra as
most of the CIB and CO emission comes from halos with masses between
$10^{12}-10^{13}M_{\odot}$ Maniyar_2021 . We then proceed to calculate the
auto- and cross-power spectra of the CIB, tSZ, and CO lines. This is shown in
Fig. 2 where we plot the power spectrum of all the signals at 150 and 220 GHz
in the left and right panels respectively for SPT. Please note that the
‘Radio’ and ‘kSZ’ power spectra which we do not simulate here are taken from
the SPT results from Reichardt_2021 . As expected, the CIB anisotropy power
spectrum (shown in blue and including the shot noise) dominates over other
signals on such small scales with an amplitude of $\mathcal{D}_{3000}^{\rm
CIB}\sim 17\>{\rm\mu K}^{2}$ at 150 GHz. The next dominant component at 150
GHz is the tSZ power spectrum (violet) with $\mathcal{D}_{3000}^{\rm tSZ}\sim
6.15\>{\rm\mu K}^{2}$. Since both the CIB and tSZ are excellent tracers of the
large-scale structure (LSS) of the Universe, they are expected to be
correlated Maniyar_2021 . This is exactly what we find here with the CIB-tSZ
cross-correlation (cyan) amplitude $\mathcal{D}_{3000}^{\rm tSZ\times
CIB}\sim-2.10\>{\rm\mu K}^{2}$ at 150 GHz. We plot the absolute value of the
CIB-tSZ cross-correlation since it is negative at 150 GHz. It is worth noting
that the SPT collaboration reported the corresponding values for the total
CIB, tSZ, and CIB-tSZ power spectra to be $\mathcal{D}_{3000}\sim
12,3.42,-0.9\>{\rm\mu K}^{2}$ respectively for 150 GHz. At 220 GHz, we get a
CIB power spectrum amplitude of $115\>{\rm\mu K}^{2}$ where SPT collaboration
find a value of $121.3\>{\rm\mu K}^{2}$. As expected, since the 220 GHz
channel is at the center of tSZ null, both the tSZ and its cross-power
spectrum with the CIB is negligible. Our results are therefore comparable with
the SPT results.

Although we present our results for 150 and 220 GHz SPT channels here, similar
analysis can be carried out at 95 GHz channel as well as for the cross-power
spectra between all these frequencies.

###  III.2 CO lines and their cross-correlations

Here we discuss the power spectrum of the CO lines and their cross-
correlations with other signals. Instead of showing the individual
CO($J\rightarrow J-1$) lines, in Fig. 2, we plot the sum of their power
spectra for clarity. As can be seen, at $\ell=3000$ the sum of all the CO line
power spectra, shown in blue, is $\mathcal{D}_{3000}^{\rm CO,total}\sim
0.08\>{\rm\mu K}^{2}$ at 150 GHz. In comparison to the other signals, CO power
spectrum is at least an order of magnitude smaller at 150 GHz. For 220 GHz, we
find that the total CO power spectrum has the same amplitude as that of the
radio component. However, at $\ell=3000$, they are both much smaller than the
CIB and kSZ power. This therefore supports the previous CMB power spectrum
analysis studies which did not account for the presence of the CO lines in the
CMB maps.

However, things get interesting when we consider the cross-correlations
between CO and the CIB and tSZ. CO lines are tracers of star-forming galaxies
(and hence the LSS) similar to the CIB. They are therefore expected to be
correlated with the CIB and tSZ. This is exactly what we find. For the CIB-CO
cross-power spectrum, we get $\mathcal{D}_{3000}^{\rm CIB\times CO,total}\sim
0.89\>{\rm\mu K}^{2}$ at 150 GHz. The CIB-CO cross-correlation has therefore
similar power as (absolute) CIB-tSZ cross-correlation at 150 GHz. On the other
hand, we get $\mathcal{D}_{3000}^{\rm CIB\times CO,total}\sim 3.20\>{\rm\mu
K}^{2}$ at 220 GHz making it the biggest signal after the CIB. This shows that
although we can neglect the auto-power spectrum of different CO lines due to
their very small amplitude, their cross-correlation with the CIB can be large
and significant already for the current data.

Figure 3: Two-dimension contour plot showing the fraction ($10\%$ to $90\%$)
of the total tSZ power at $\ell=3000$ as a function of different mass and
redshift ranges. This shows that the majority ($\sim 70\%$) of the tSZ power
is sourced by massive halos $M_{200}>10^{13}M_{\odot}$ at low redshifts $z<1$.

On the other hand, although the tSZ signal is quite large on these scales, we
find that the cross-correlation between the tSZ and CO lines is small and not
as important as the CIB-CO correlation. This is expected at 220 GHz where the
tSZ signal is itself negligible. At 150 GHz, its amplitude
$\mathcal{D}_{3000}^{\rm tSZ\times CO,total}\sim-0.18\>{\rm\mu K}^{2}$ is
similar to the total CO power spectrum. In order to understand this, we
explore the origin of the tSZ power. Similar to Bhattacharya_2012 , we plot
the fractional tSZ power for halos with mass $>M_{200}$ and redshift $<z$,
from 10-90% of the total power spectrum from our simulation in Fig. 3. This
plot shows that majority of the tSZ power ($\sim 70\%$) comes from massive
halos ($M_{200}>10^{13}M_{\odot}$) at low redshifts ($z<1$).

As we precisely know the rest frame frequency of different CO($J\rightarrow
J-1$) transitions, we can determine their redshift distribution within the SPT
bandpass. Since we have created our CO maps assuming a spectral bandwidth of
$\Delta\nu=50$ centered at 150 GHz, we have different CO lines at various
redshifts which fall within the observed frequencies of $130<\nu<170$ GHz.
This means that CO($1\rightarrow 0$) with the rest frame frequency of 115.27
GHz, does not show up in our maps. In fact, apart from CO($2\rightarrow 1$)
transitions from $0.35<z<0.75$, the rest of the CO lines in our map are from
higher redshifts ($z>1$). In Fig. 4, we show the power spectra of individual
CO lines from our maps. As we can see, at $\ell=3000$, the CO($5\rightarrow
4$) line ($2.4<z<3.4$) has the highest power followed by CO($4\rightarrow 3$)
line ($1.7<z<2.5$), CO($2\rightarrow 1$) line ($0.35<z<0.75$), and
CO($3\rightarrow 2$) line ($1.0<z<1.7$) which have almost the same power.
Therefore, most of the CO power comes from redshifts $z>1$. All these factors
help us understand the smaller amplitude of the tSZ-CO cross-correlation
compared to the CIB-CO cross-correlation for the SPT 150 GHz channel.

Figure 4: Power spectra of individual CO lines observed in the SPT 150 GHz
maps. Most of the power comes from the different CO lines between $1<z<3.5$.

In Appendix A, we show the cross-correlation of individual CO lines with the
CIB and tSZ in Fig. 5 and Fig. 6 respectively and discuss their implications.

##  IV Implications of this new foreground

###  IV.1 Implications for kSZ measurements

Measuring the kSZ power spectrum from the data is an important goal for the
current and upcoming CMB experiments. The SPT collaboration has performed
multi-component analyses of the small scale CMB power spectrum accounting for
different foregrounds George_2015 ; Reichardt_2021 . In fact, their recent
analysis reported a first measurement of the kSZ power at $\geq 3\sigma$ with
a value of $\mathcal{D}_{3000}^{\rm kSZ}=3.0\pm 1.0\>{\rm\mu K}^{2}$ (Fig. 2).
They also show that the kSZ power is degenerate with the corresponding tSZ and
CIB-tSZ cross-correlation power respectively. In this analysis, they do not
account for the presence of CO lines in the SPT maps. As shown in Sect. III.2,
the total auto-power spectrum of the CO lines is an order of magnitude smaller
($\sim 0.1\>{\rm\mu K}^{2}$) than the expected kSZ power spectrum ($\sim
1\>{\rm\mu K}^{2}$), so this should not affect the results. However, here we
have shown that the cross-correlation of the CIB with CO lines has an
amplitude similar to the kSZ and CIB-tSZ cross-correlation. Looking only at
the 150 GHz, a significant (negative) cross-correlation between the tSZ and CO
lines could have resulted in a fortuitous cancellation of the CIB-CO
correlation, which is not the case here. This is a very important result and
motivates current and future analyses to include this CIB-CO correlation
component.

Since this new signal has almost the same power as the kSZ, it may
significantly bring down the measured kSZ power. On the other hand, the shape
of this CIB-CO correlation seems similar to the other signals such as the CIB
and CIB-tSZ correlation. This is different than the kSZ which is mostly flat
on these scales. It is thus possible that without affecting the kSZ power,
this new signal may change values of the CIB and CIB-tSZ correlation (but also
other components as radio). Without performing a thorough multi-component
analysis, it is therefore not possible to say whether this will decrease the
amount or signal-to-noise ratio of the kSZ power detected by the SPT team. In
an upcoming paper, we will perform such an analysis to determine the impact of
CIB-CO correlations on kSZ power spectrum.

###  IV.2 CO as a signal

CO($J\rightarrow J-1$) molecular line emissions are among the brightest lines
in the galaxy spectra. Their measurements can be used to estimate the amount
of gas mass and star formation rate in galaxies Bernal_2022 . They also act as
a proxy for cosmic molecular gas history. Thus CO lines are one of the prime
targets of some of the current and upcoming LIM surveys. Recently,
Keating_2016  and Keating_2020  reported a $2\sigma$ and $4\sigma$ detection
of the shot-noise component of the CO power spectrum. However, to date, there
has been no detection of the clustering component of the CO power spectrum. As
we showed in the previous section, the CIB-CO correlation acts as a
significant signal in the small-scale CMB data. Thus, those data can be
potentially used to measure the CIB-CO cross-correlation signal for the first
time.

As mentioned in Sect. I, predictions for CO power spectra can vary by a large
amount between different models. Such a measurement will therefore not only
help in ruling out some existing models but also constrain the parameter space
of the most promising models. This will, in turn, enable us to draw
conclusions on e.g. the cosmic molecular gas history. As several CO lines from
different redshifts show up at different frequencies in the CMB maps,
measuring this amplitude across these frequency channels has the potential to
break some potential degeneracies in the models.

###  IV.3 Implications for studies involving cross-correlations

In view of the high quality data from the current and upcoming cosmological
experiments, cross-correlations of different tracers of the LSS have emerged
as a powerful tool to get rid of different biases and isolate signals of
interest. All such analyses involving small-scale CMB information will be
affected by this new signal. One example is the so-called projected-field
estimator to detect the kSZ signal Ferraro_2018 . It involves squaring the CMB
map and cross-correlating it with a LSS tracer or in general measuring the
$\langle TT\delta\rangle$ three-point function where $T$ is the CMB data and
$\delta$ is a LSS tracer. Such analyses will be affected by the CO lines
present in the CMB maps as they will cross-correlate with the LSS tracer and
CIB, and therefore act as a source of bias (as $\langle{\rm CO}-{\rm
CO}-\delta\rangle$ and $\langle{\rm CO}-{\rm CIB}-\delta\rangle$) terms to the
kSZ signal if not properly accounted for. In fact, this could be one of the
components responsible for an unexplained excess kSZ amplitude observed by
Kusiack_2021  using this estimator. Therefore, with the more precise CMB data,
CO emission in galaxies may give non-significant contributions and thus must
be accounted for in cross-correlation analyses.

##  V Conclusion

As high-resolution CMB data keeps on getting more precise with the current and
upcoming experiments, this presents us with an exciting opportunity to perform
a wide range of analyses. In order to best utilize this rich data, it is
therefore crucial to understand and account for all the sources contributing
at these millimeter frequencies. In this paper, we study the extragalactic CO
molecular emission lines in such a context, which traditionally have been
neglected in past analyses.

Since the CMB experiments use broadband filters, we thus expect different CO
lines from several redshifts to fall within these bands. Using the SIDES
simulations, we prepare a map containing the CIB, tSZ, and CO lines as
observed by the SPT telescope at 150 and 220 GHz. These are the first such
simulations consistently containing all these components. We show that our
predictions of the CIB, tSZ power spectra, and their cross-correlations match
the observations. We show that the power spectrum of the CO lines is an order
of magnitude smaller than the other foregrounds. Their contribution to the
total power spectrum at these scales is thus negligible. However, we then go
on to show that the cross-correlation of the CIB with CO lines has a
contribution similar to the CIB-tSZ correlation and the kSZ power spectrum. It
therefore contributes a non-negligible amount to the total power at these
scales. This is the main result of this paper. We also show that, unlike the
CIB-CO, the tSZ-CO correlation is not as strong at 150 GHz. This is due to the
difference in the redshift and halo masses of the sources contributing to
these two signals. Similar arguments are expected to hold at 95 GHz as well
which is the other main CMB channel.

Our result is of special significance in light of the recently reported first-
ever $>3\sigma$ detection of the kSZ power spectrum by Reichardt_2021 . This
result is quite important as kSZ measurements and their interrogations promise
to shed light on the reionization of the Universe. They perform this
measurement of the kSZ accounting for the presence of several other
foregrounds like the CIB, tSZ, and radio emission. However, they do not
account for the CO lines and their cross-correlations with the CIB which we
show here have an amplitude similar to the kSZ power spectrum. As discussed in
Sect. IV.1, our results call for a more careful analysis of the data
accounting for CIB-CO cross-correlations, which may have a significant impact
on the kSZ measurements.

Our results also show that the CO lines present in the CMB maps will be
crucial for all the small scale CMB auto-power spectrum and cross-correlation
studies involving a LSS tracer. Upcoming CMB experiments aim to use small
scale CMB measurements to constrain a diverse range of physics including
precision constraints on the $\Lambda$CDM model to neutrino mass to decaying
dark matter (e.g. Abazajian_16, ). The signatures of many of these processes
are very subtle and a preeminent challenge is mitigating possible biases. Not
modelling the CO lines risks their power being incorrectly attributed to novel
physics. Hence, including the CO lines in these analysis is essential to avoid
any potential biases.

In an upcoming paper, we will perform a comprehensive study of the small-scale
CMB power spectrum with all these foregrounds. This analysis will not only
determine the impact of CIB-CO correlations on kSZ power spectrum measurements
but will also provide us with the first CIB-CO cross-correlation measurement
from the CMB maps. This measurement will help us constrain the star formation
history of the Universe as both the CIB and CO trace star-forming galaxies
within our Universe.

Finally, in this paper, we solely focused on the extragalactic CO molecular
lines. There are other lines like [CI] contributing to the millimeter sky
Bethermin_2022 . As the data quality keeps on getting better, we might have to
look beyond CO lines and also account for other lines like [CI].

###### Acknowledgements.

We thank Srinivasan Raghunathan from the SPT collaboration for providing us
with the SPT measurements of the Radio and kSZ power spectrum. We thank
Matthieu Béthermin for his help in using SIDES. This project has received
funding from the European Research Council (ERC) under the European Union’s
Horizon 2020 research and innovation program (grant agreement No 788212) and
from the Excellence Initiative of Aix-Marseille University-A*Midex, a French
“Investissements d’Avenir” program. Some of the results in this paper have
been derived using various standard Python packages (numpy, scipy,
matplotlib). Part of this work was carried out on the ancestral land of the
Muwekma Ohlone Tribe, successors of the Verona Band of Alameda County.

## References

  * (1) Planck Collaboration et al., A&A641, A6 (2020).
  * (2) Planck Collaboration et al., A&A594, A22 (2016).
  * (3) J. Dunkley et al., J. Cosmology Astropart. Phys2013, 025 (2013).
  * (4) E. M. George et al., ApJ799, 177 (2015).
  * (5) P. Ade et al., J. Cosmology Astropart. Phys2019, 056 (2019).
  * (6) K. N. Abazajian et al., arXiv e-prints (2016).
  * (7) C. L. Reichardt et al., ApJ908, 199 (2021).
  * (8) E. D. Kovetz et al., arXiv e-prints arXiv:1709.09066 (2017).
  * (9) J. L. Bernal and E. D. Kovetz, arXiv e-prints arXiv:2206.15377 (2022).
  * (10) M. Righi, C. Hernández-Monteagudo, and R. A. Sunyaev, A&A489, 489 (2008).
  * (11) H. Padmanabhan, MNRAS475, 1477 (2018).
  * (12) G. K. Keating et al., ApJ830, 34 (2016).
  * (13) M. Béthermin et al., A&A 607, A89 (2017).
  * (14) M. Béthermin et al., A&A667, A156 (2022).
  * (15) A. Gkogkou et al., arXiv e-prints arXiv:2212.02235 (2022).
  * (16) T. Ishiyama et al., Monthly Notices of the Royal Astronomical Society 506, 4210 (2021).
  * (17) I. Davidzon et al., ApJ 852, 107 (2018).
  * (18) C. Schreiber et al., A&A 575, A74 (2015).
  * (19) R. C. Kennicutt, Annual Review of Astronomy and Astrophysics 36, 189 (1998).
  * (20) G. E. Magdis et al., ApJ 760, 6 (2012).
  * (21) M. Béthermin et al., A&A 573, A113 (2015).
  * (22) M. T. Sargent et al., ApJ 793, 19 (2014).
  * (23) F. Bournaud et al., A&A 575, A56 (2015).
  * (24) E. Daddi et al., A&A 577, A46 (2015).
  * (25) G. Stein et al., J. Cosmology Astropart. Phys2020, 012 (2020).
  * (26) N. Sehgal et al., ApJ709, 920 (2010).
  * (27) H. Trac, P. Bode, and J. P. Ostriker, ApJ727, 94 (2011).
  * (28) K. Osato and D. Nagai, arXiv e-prints arXiv:2201.02632 (2022).
  * (29) N. Battaglia, J. R. Bond, C. Pfrommer, and J. L. Sievers, ApJ758, 75 (2012).
  * (30) G. Lagache et al., A&A642, A232 (2020).
  * (31) B. Bolliet, B. Comis, E. Komatsu, and J. F. Macías-Pérez, MNRAS477, 4957 (2018).
  * (32) S. K. Choi et al., J. Cosmology Astropart. Phys2020, 045 (2020).
  * (33) E. Komatsu and U. Seljak, MNRAS336, 1256 (2002).
  * (34) G. Stein et al., J. Cosmology Astropart. Phys2020, 012 (2020).
  * (35) A. Maniyar, M. Béthermin, and G. Lagache, A&A645, A40 (2021).
  * (36) S. Bhattacharya et al., ApJ760, 5 (2012).
  * (37) G. K. Keating, D. P. Marrone, G. C. Bower, and R. P. Keenan, ApJ901, 141 (2020).
  * (38) S. Ferraro and J. C. Hill, Phys. Rev. D97, 023512 (2018).
  * (39) A. Kusiak et al., Phys. Rev. D104, 043518 (2021).
  * (40) A. S. Maniyar, M. Béthermin, and G. Lagache, A&A614, A39 (2018).

##  Appendix A Correlations of individual CO lines with foregrounds

Fig. 5 shows the cross-power spectra between different CO lines and the CIB at
150 GHz. The correlation is highest for CO($5\rightarrow 4$) line from
$2.4<z<3.4$. This is followed by CO($4\rightarrow 3$) and, CO($3\rightarrow
2$) lines. This is similar to the individual CO line power spectra plotted in
Fig. 4. However, in Fig. 4, CO($2\rightarrow 1$) has the same amplitude as
that of CO($3\rightarrow 2$) which is not the case here where CIB-
CO($3\rightarrow 2$) correlation is much higher than CIB-CO($2\rightarrow 1$).
This can be explained as follows. As the CIB is traced by star-forming
galaxies, bulk of the CIB emission comes from the epoch when the star
formation in the Universe peaked between $1<z<3$ Maniyar_2018 . Also, the CIB
has this peculiar property that low frequency CIB maps trace galaxies at
higher redshifts and vice-versa. This also explains higher correlation of the
CIB with the CO($3\rightarrow 2$) line than CO($2\rightarrow 1$).

Figure 5: Cross-power spectra of individual CO lines with the CIB at 150 GHz.

On the other hand, as can be seen from Fig. 6, this is not the case for tSZ-CO
correlations. Although CO($5\rightarrow 4$) line has the highest auto-power
spectrum at 150 GHz, tSZ has the highest correlation with CO($2\rightarrow 1$)
line originating between $0.4<z<0.8$. As explained in Sect. III.2, this is due
to the low redshift origin of most of the tSZ power. In fact, we can see from
Fig. 3 that most of the tSZ power comes from massive halos at low redshifts
($z<1$). Therefore, as the higher CO transition lines trace higher redshifts,
their correlation with the tSZ keeps going down as well.

Figure 6: Cross-power spectra of individual CO lines with the tSZ at 150 GHz.


# On the Interpolating Sesqui-Harmonicity of Vector Fields

Bouazza Kacimi  University of Mascara
Faculty of Exact Sciences
Department of Mathematics
BP 305 Route de Mamounia
29000 Mascara
Algeria<EMAIL_ADDRESS>Amina Alem  University of Mascara
Faculty of Exact Sciences
Department of Mathematics
BP 305 Route de Mamounia
29000 Mascara
Algeria<EMAIL_ADDRESS>Mustafa Özkan  Gazi University
Faculty of Sciences
Department of Mathematics
06500 Ankara
Turkey<EMAIL_ADDRESS>
(Date: January 1, 2004)

###### Abstract.

This article deals with the interpolating sesqui-harmonicity of a vector field
$X$ viewed as a map from a Riemannian manifold $(M,g)$ to its tangent bundle
$TM$ endowed with the Sasaki metric $g_{S}$. We show characterization theorem
for $X$ to be interpolating sesqui-harmonic map. We give also the critical
point condition which characterizes interpolating sesqui-harmonic vector
fields. When $(M,g)$ is compact and oriented and under some conditions, we
prove that $X$ is an interpolating sesqui-harmonic vector field (resp.
interpolating sesqui-harmonic map) if and only if $X$ is parallel. Moreover,
we extend this result for a left-invariant vector field on a Lie group $G$
having a discrete subgroup $\Gamma$ such that the quotient $\Gamma\backslash
G$ is compact.

###### Key words and phrases:

Tangent bundle, Sasaki metric, Interpolating sesqui-harmonic maps

###### 1991 Mathematics Subject Classification:

Primary 58E20; Secondary 53C20

##  1\. Introduction

The theory of harmonic maps is of essential importance in both geometry,
analysis and physics. For example, on mathematical side harmonic maps are
among the most studied variational problems in geometric analysis, and on
physical side they are closely related to nonlinear field theory in
theoretical physics. Given a smooth map $\varphi:M\rightarrow N$ between two
Riemannian manifolds $(M,g)$ and $(N,h)$ of dimensions $m$ and $n$
respectively, the map $\varphi$ is said to be harmonic if it is a critical
point of the energy functional defined by:

$E(\varphi)=\frac{1}{2}\int_{D}\left\|d\varphi\right\|^{2}v_{g},$ (1.1)

for any compact domain $D$ of $M$, where $v_{g}$ is the volume element of
$(M,g)$. The Euler-Lagrange equation of $E(\varphi)$ is [2, 9]

$\tau(\varphi)=\operatorname{Tr}_{g}(\nabla
d\varphi)=\sum_{i=1}^{m}\\{\nabla^{\varphi}_{e_{i}}d\varphi(e_{i})-d\varphi(\nabla_{e_{i}}e_{i})\\}=0,$

where $\tau(\varphi)$ is the tension field of $\varphi$, here
$\nabla^{\varphi}$ is the connection on the vector bundle $\varphi^{-1}TN$
induced from the Levi-Civita connection $\nabla^{N}$ of $(N,h)$, $\nabla$ is
the Levi-Civita connection of $(M,g)$ and $\\{e_{i}\\}_{i=1}^{m}$ is a local
orthonormal frame field of $(M,g)$. For a recent survey on harmonic maps see
[12]. Among the first generalizations of harmonic maps is the notion of
polyharmonic maps of order $k$ between Riemannian manifolds introduced by
Eells and Lemaire in [8]. For $k=2$, they defined the bienergy of $\varphi$ as
the functional

$E_{2}(\varphi)=\frac{1}{2}\int_{D}\left\|\tau(\varphi)\right\|^{2}v_{g},$
(1.2)

for any compact domain $D$ of $M$. The map $\varphi$ is said to be biharmonic
if it is a critical point of the bienergy functional (1.2). The associated
Euler-Lagrange equation is derived by Jiang [14] as follows:

$\tau_{2}(\varphi)=\Delta^{\varphi}\tau(\varphi)-\overset{m}{\underset{i=1}{\sum}}R^{N}(\tau(\varphi),d\varphi(e_{i}))d\varphi(e_{i})=0,$

where $\tau_{2}(\varphi)$ is the bitension field of $\varphi$,
$\Delta^{\varphi}\tau(\varphi)=-\sum_{i=1}^{m}(\nabla^{\varphi}_{e_{i}}\nabla^{\varphi}_{e_{i}}\tau(\varphi)-\nabla^{\varphi}_{\nabla^{M}_{e_{i}}e_{i}}\tau(\varphi))$
is the rough Laplacian on $\varphi^{-1}TN$, and $R^{N}$ is the curvature
tensor of the target manifold $N$. By definition, it can be seen that every
harmonic map is biharmonic. However, a biharmonic map can be non-harmonic in
which case it is called proper biharmonic. Biharmonic maps have been studied
by several researchers see for instance [20, 22] and references therein.

Branding [4], introduced an action functional for maps between Riemannian
manifolds that interpolated between the actions for harmonic and biharmonic
maps as follows:

$E_{\delta_{1},\delta_{2}}(\varphi)=\delta_{1}\int_{D}\left\|d\varphi\right\|^{2}v_{g}+\delta_{2}\int_{D}\left\|\tau(\varphi)\right\|^{2}v_{g},$
(1.3)

where $\delta_{1},\delta_{2}\in\mathbb{R}$. The functional (1.3) is used in
string theory of physics, it is known as bosonic string with extrinsic
curvature term, see [21]. The map $\varphi$ is said to be interpolating
sesqui-harmonic if it is a critical point of the functional (1.3). The
associated Euler-Lagrange equation is [4]

$\tau_{\delta_{1},\delta_{2}}(\varphi)=\delta_{1}\tau(\varphi)+\delta_{2}\tau_{2}(\varphi)=0.$
(1.4)

We mention that the field $\tau_{\delta_{1},\delta_{2}}(\varphi)$ is appeared
with different sign from the expression in [4] due to the different sign of
the bitension field. It is obvious that harmonic maps solve (1.4). Moreover,
if $\delta_{1}=0$ and $\delta_{2}=1$, then an interpolating sesqui-harmonic
map turns into a biharmonic map. So, harmonic maps and biharmonic maps are
subclass of interpolating sesqui-harmonic maps. Note that there have been
several articles dealing with particular aspects of (1.3) see [4]. In [5],
Branding studied various analytic aspects of interpolating sesqui-harmonic
maps between Riemannian manifolds particulary a spherical target. In [15],
Karaca et al. considered interpolating sesqui-harmonic curves in Sasakian
space forms and obtained the necessary and sufficient conditions for Legendre
curves to be interpolating sesqui-harmonic.

On the other hand, denote by $\mathfrak{X}(M)$ the set of all smooth vector
fields on $M$ and by $g_{S}$ the Sasaki metric on the tangent bundle $TM$. Any
$X\in\mathfrak{X}(M)$ defines a smooth map from $(M,g)$ to $(TM,g_{S})$. When
$M$ is compact, it was proved in [13] and also in [19] that
$X:(M,g)\rightarrow(TM,g_{S})$ is an harmonic map if and only if $X$ is
parallel; moreover this result remains true if $X$ is a harmonic section of
$TM$, i.e., a critical point of the energy functional $E$ restricted to the
set $\mathfrak{X}(M)$ [10]. The bienergy of $X\in\mathfrak{X}(M)$ is the
bienergy of the corresponding map see [17]. In [17], the authors obtained the
critical point of the bienergy functional $E_{2}$ restricted to the set
$\mathfrak{X}(M)$ (equivalently, $X$ is a biharmonic vector field) and showed
that if $M$ is compact, then $X$ is biharmonic vector field (resp. biharmonic
map) if and only if $X$ is parallel. In [1], the authors established the
formula of the bitension field of $X:(M,g)\rightarrow(TM,g_{S})$, and obtained
characterization theorem for $X$ to be biharmonic map, furthermore they gave
the relationship between the notion of biharmonic vector field and that of
vector field which is biharmonic map. Especially, they proved that a left-
invariant vector field $X$ on three dimensional unimodular Lie group is
biharmonic vector field (resp. biharmonic map) if and only if $X$ is parallel.

In this note, we will study the interpolating sesqui-harmonicity of
$X\in\mathfrak{X}(M)$ viewed as a map $X:(M,g)\rightarrow(TM,g_{S})$ building
on the developed results for harmonic vector fields [7, 10], vector fields
which are harmonic maps [13, 19], biharmonic vector fields [17] and vector
fields which are biharmonic maps [1]. So, it is natural to address the problem
of characterizing those vector fields for which the corresponding map is an
interpolating sesqui-harmonic map, and also to check whether vector fields $X$
that are critical points of the functional $E_{\delta_{1},\delta_{2}}$
restricted to variations through vector fields, such a vector field is called
interpolating sesqui-harmonic vector field.

This paper is organized as follows. Section $2$ contains some basic notions,
properties and results that will be needed later. In section $3$, we establish
the formula of the field $\tau_{\delta_{1},\delta_{2}}(X)$ of
$X:(M,g)\rightarrow(TM,g_{S})$ (see Theorem 3.1) and provide the conditions
which characterize vector fields which are interpolating sesqui-harmonic maps
(see Theorem 3.2). By means of the formula of the field
$\tau_{\delta_{1},\delta_{2}}(X)$ of $X$, we derive the first variational
formula associated to the energy functional $E_{\delta_{1},\delta_{2}}$
restricted to the space $\mathfrak{X}(M)$ (see Theorem 3.4). As a corollary,
we obtain the critical point condition characterizes interpolating sesqui-
harmonic vector field (see Corollary 3.6). Consequently, we obtain the
relationship between interpolating sesqui-harmonic vector fields and vector
fields which are interpolating sesqui-harmonic maps (see Corollary 3.7).
Afterwards, if $\delta_{1}$ and $\delta_{2}$ have the same sign, and $(M,g)$
is compact and oriented, we prove that $X$ is an interpolating sesqui-harmonic
vector field (resp. interpolating sesqui-harmonic map) if and only if $X$ is
parallel (see Theorem 3.8 and Theorem 3.9), we present an example of non-
parallel vector field which is an interpolating sesqui-harmonic map on the Sol
space (see Example 3.10). Finally, in section $4$, if $\delta_{1}$ and
$\delta_{2}$ have the same sign, and $G$ is a Lie group having a discrete
subgroup $\Gamma$ such that the quotient $\Gamma\backslash G$ is compact, we
prove that a left-invariant vector field $X$ on $G$ is interpolating sesqui-
harmonic (resp. interpolating sesqui-harmonic map) if and only if $X$ is
parallel (see Theorem 4.1). Thereby, when $\delta_{1}$ and $\delta_{2}$ have
different signs, we completely determine the set of left-invariant
interpolating sesqui-harmonic vector fields and left-invariant vector fields
which are interpolating sesqui-harmonic maps on the Heisenberg group $Nil$
(see Theorem 4.6).

##  2\. Preliminaries

We recall here some basic facts on the geometry of tangent bundle. We refer
the reader to [6, 16, 23] and references therein for further details. Let
$(M,g)$ be a Riemannian manifold of dimension $m$ and $(TM,\pi,M)$ be its
tangent bundle, a local chart $(U,x^{i})_{1\leq i\leq m}$ on $M$ induces a
local chart $(\pi^{-1}(U),x^{i},y^{i})_{1\leq i\leq m}$ on $TM$. Denotes the
Christoffel symbols of $g$ by $\Gamma_{jk}^{i}$, the tangent space
$T_{(x,u)}TM$ at a point $(x,u)$ in $TM$ is a direct sum of the vertical
subspace $\mathcal{V}_{(x,u)}=\ker(d\pi\mid_{(x,u)})$ and the horizontal
subspace $\mathcal{H}_{(x,u)}$, with respect to the Levi-Civita connection
$\nabla$ of $M$:

$T_{(x,u)}TM=\mathcal{H}_{(x,u)}\oplus\mathcal{V}_{(x,u)}.$

Let $X|_{U}=X^{i}\frac{\partial}{\partial x^{i}}$ be a local vector field on
$M$. The vertical and the horizontal lifts of $X$ are defined respectively by:

$X^{v}|_{\pi^{-1}(U)}=(X^{i}\circ\pi)\frac{\partial}{\partial y^{i}},$

and

$X^{h}|_{\pi^{-1}(U)}=(X^{i}\circ\pi)\frac{\partial}{\partial
x^{i}}-(\Gamma_{jk}^{i}\circ\pi)(X^{j}\circ\pi)y^{k}\frac{\partial}{\partial
y^{i}}.$

The Sasaki metric on $TM$ is the Riemannian metric $g_{S}$ defined by

$g_{S}(X^{h},Y^{h})=g_{S}(X^{v},Y^{v})=g(X,Y)\circ\pi,\;g_{S}(X^{v},Y^{h})=0,$
(2.1)

for any $X,Y\in\Gamma(TM)$.

A vector field $X$ on $(M,g)$ can be viewed as the immersion
$X:(M,g)\rightarrow(TM,g_{S}):x\mapsto(x,X_{x})\in TM$ into its tangent bundle
$TM$ equipped with the Sasaki metric $g_{S}$. If $Y\in\Gamma(TM)$ then, we
have (see [7, pp. 50])

$dX(Y)=\\{Y^{h}+(\nabla_{Y}X)^{v}\\}\circ X.$ (2.2)

The tension field $\tau(X)$ is given by [10] and rewritten in [17] as follows:

$\tau(X)=(-S(X))^{h}+(-\bar{\Delta}X)^{v},$ (2.3)

where $S(X)=\overset{m}{\underset{i=1}{\sum}}R(\nabla_{e_{i}}X,X)e_{i}$, here
$R$ denotes the Riemannian curvature tensor taken with the sign convention

$R(X,Y)Z=\nabla_{X}\nabla_{Y}Z-\nabla_{Y}\nabla_{X}Z-\nabla_{[X,Y]}Z,$

for all vector fields $X$, $Y$ and $Z$ on $M$, and $\bar{\Delta}X$ is the
rough Laplacian given by
$\bar{\Delta}X=-tr_{g}(\nabla^{2}X)=\overset{m}{\underset{i=1}{\sum}}(\nabla_{\nabla_{e_{i}}e_{i}}X-\nabla_{e_{i}}\nabla_{e_{i}}X).$
A vector field $X$ defines a harmonic map from $(M,g)$ to $(TM,g_{S})$ if and
only if $\tau(X)=0$, equivalently $\bar{\Delta}X=0$ and $S(X)=0$, $X$ is
called harmonic vector field if it is a critical point of the energy
functional (1.1), only considering variations through vector fields. The
corresponding Euler-Lagrange equation is given by $\bar{\Delta}X=0,$ so $X$ is
a harmonic map if and only if $X$ is a harmonic vector field and $S(X)=0$.

##  3\. The interpolating sesqui-harmonicity of vector fields

In the sequel, we give the formula of the field
$\tau_{\delta_{1},\delta_{2}}(X)$ of $X:(M,g)\rightarrow(TM,g_{S})$.

######  Theorem 3.1.

Let $(M,g)$ be an $m$-dimensional Riemannian manifold and $(TM,g_{S})$ its
tangent bundle equipped with the Sasaki metric, if
$X:(M,g)\rightarrow(TM,g_{S})$ is a smooth vector field then the field
$\tau_{\delta_{1},\delta_{2}}(X)$ of $X$ is given by

$\displaystyle\tau_{\delta_{1},\delta_{2}}$
$\displaystyle(X)=\Big{\\{}-\delta_{1}\bar{\Delta}X-\delta_{2}\bar{\Delta}\bar{\Delta}X-\delta_{2}\sum_{i=1}^{m}[(\nabla_{e_{i}}R)(e_{i},S(X))X$
$\displaystyle+R(e_{i},\nabla_{e_{i}}S(X))X+2R(e_{i},S(X))\nabla_{e_{i}}X]\Big{\\}}^{v}+\Big{\\{}-\delta_{1}S(X)-\delta_{2}\bar{\Delta}S(X)$
$\displaystyle-\delta_{2}R(X,\bar{\Delta}X)S(X)+\delta_{2}\sum_{i=1}^{m}[R(X,\nabla_{e_{i}}\bar{\Delta}X)e_{i}-R(\nabla_{e_{i}}X,\bar{\Delta}X)e_{i}$
$\displaystyle-R(e_{i},S(X))e_{i}-(\nabla_{S(X)}R)(\nabla_{e_{i}}X,X)e_{i}+R(X,\nabla_{e_{i}}X)\nabla_{e_{i}}S(X)$
$\displaystyle-R(X,R(e_{i},S(X))X)e_{i}]\Big{\\}}^{h}.$ (3.1)

###### Proof.

On making use of (1.4) and (2.3) and Theorem 3.1 in [1], we find easily the
formula (3.1). ∎

Then, we deduce

######  Theorem 3.2.

Let $(M,g)$ be an $m$-dimensional Riemannian manifold and
$X\in\mathfrak{X}(M)$, then $X:(M,g)\rightarrow(TM,g_{S})$ is an interpolating
sesqui-harmonic map if and only if

$\delta_{1}\bar{\Delta}X+\delta_{2}\bar{\Delta}\bar{\Delta}X+\delta_{2}\sum_{i=1}^{m}[(\nabla_{e_{i}}R)(e_{i},S(X))X\\\
+R(e_{i},\nabla_{e_{i}}S(X))X+2R(e_{i},S(X))\nabla_{e_{i}}X]=0,$

and

$\delta_{1}S(X)+\delta_{2}\bar{\Delta}S(X)+\delta_{2}R(X,\bar{\Delta}X)S(X)-\delta_{2}\sum_{i=1}^{m}[R(X,\nabla_{e_{i}}\bar{\Delta}X)e_{i}\\\
-R(\nabla_{e_{i}}X,\bar{\Delta}X)e_{i}-R(e_{i},S(X))e_{i}-(\nabla_{S(X)}R)(\nabla_{e_{i}}X,X)e_{i}\\\
+R(X,\nabla_{e_{i}}X)\nabla_{e_{i}}S(X)-R(X,R(e_{i},S(X))X)e_{i}]=0,$

where $\\{e_{i}\\}_{i=1}^{m}$ is a local orthonormal frame field of $(M,g)$.

Let $(M,g)$ be a compact $m$-dimensional Riemannian manifold. The energy
functional $E_{\delta_{1},\delta_{2}}$ of $X$ is defined to be the energy of
the corresponding map $X:(M,g)\longrightarrow(TM,g_{S})$. More precisely,
combining relations (2.1), (2.2) and (2.3), one obtain

$\displaystyle E_{\delta_{1},\delta_{2}}(X)$
$\displaystyle=\delta_{1}\int_{M}\left\|dX\right\|^{2}v_{g}+\delta_{2}\int_{M}\left\|\tau(X)\right\|^{2}v_{g},$
$\displaystyle=\delta_{1}m\;Vol(M)+\delta_{1}\int_{M}\left\|\nabla
X\right\|^{2}v_{g}+\delta_{2}\int_{M}\big{[}\left\|S(X)\right\|^{2}+\left\|\bar{\Delta}X\right\|^{2}\big{]}v_{g}.$

######  Definition 3.3.

Let $(M,g)$ be a Riemannian manifold. A vector field $X\in\mathfrak{X}(M)$ is
called interpolating sesqui-harmonic if the corresponding map
$X:(M,g)\longrightarrow(TM,g_{S})$ is a critical point for the functional
$E_{\delta_{1},\delta_{2}}$, only considering variations through vector
fields.

In what follows, we determine the first variational formula of the functional
$E_{\delta_{1},\delta_{2}}$ restricted to the space $\mathfrak{X}(M)$. We
prove the following theorem:

######  Theorem 3.4.

Let $(M,g)$ be a compact oriented $m$-dimensional Riemannian manifold,
$\\{e_{i}\\}_{i=1}^{m}$ a local orthonormal frame field of $(M,g)$, $X$ a
tangent vector field on $M$ and
$E_{\delta_{1},\delta_{2}}:\mathfrak{X}(M)\longrightarrow[0,+\infty)$ the
functional $E_{\delta_{1},\delta_{2}}$ restricted to the space of all vector
fields. Then

$\displaystyle\frac{d}{dt}E_{\delta_{1},\delta_{2}}(X_{t})\bigg{|}_{t=0}=$
$\displaystyle\int_{M}\Big{\\{}g(2\delta_{1}\bar{\Delta}X+2\delta_{2}\bar{\Delta}\bar{\Delta}X+2\delta_{2}\sum_{i=1}^{m}\big{[}(\nabla_{e_{i}}R)(e_{i},S(X))X$
$\displaystyle+R(e_{i},\nabla_{e_{i}}S(X))X+2R(e_{i},S(X))\nabla_{e_{i}}X\big{]},V)\Big{\\}}v_{g},$
(3.2)

for any smooth $1$-parameter variation
$U:M\times(-\epsilon,\epsilon)\rightarrow TM$ of $X$ through vector fields
i.e., $X_{t}(z)=U(z,t)\in T_{z}M$ for any $|t|<\epsilon$ and $z\in M$, or
equivalently $X_{t}\in\mathfrak{X}(M)$ for any $|t|<\epsilon$. Also, $V$ is
the tangent vector field on $M$ given by

$V(z)=\frac{d}{dt}X_{z}(0),\quad z\in M,$

where $X_{z}(t)=U(z,t),\,(z,t)\in M\times(-\epsilon,\epsilon)$.

###### Proof.

Let $U:M\times(-\epsilon,\epsilon)\rightarrow TM$ be a smooth variation of
$X$( i.e., $U(z,0)=X(z)$ for any $z\in M$) such that $X_{t}(z)=U(z,t)\in
T_{z}M$ for any $z\in M$ and any $|t|<\epsilon$. We have

$E_{\delta_{1},\delta_{2}}(X_{t})=\delta_{1}\int_{M}\left\|dX_{t}\right\|^{2}v_{g}+\delta_{2}\int_{M}\left\|\tau(X_{t})\right\|^{2}v_{g}.$

And, from [4] one has

$\frac{d}{dt}E_{\delta_{1},\delta_{2}}(X_{t})\bigg{|}_{t=0}=-2\int_{M}g_{S}(\mathcal{V},\tau_{\delta_{1},\delta_{2}}(X))v_{g},$

where $\mathcal{V}(z)=\frac{dX_{t}(z)}{dt}\big{|}_{t=0},$ and from [7, pp.
58], one obtain

$\mathcal{V}=V^{v}\circ X.$ (3.3)

By virtue of (3.3) and the formula of $\tau_{\delta_{1},\delta_{2}}(X)$ given
by (3.1), we get

$\displaystyle\frac{d}{dt}E_{\delta_{1},\delta_{2}}(X_{t})\bigg{|}_{t=0}=$
$\displaystyle-2\int_{M}g_{S}(V^{v},\tau_{\delta_{1},\delta_{2}}(X))v_{g},$
$\displaystyle=$ $\displaystyle
2\int_{M}\Big{\\{}g(V,\delta_{1}\bar{\Delta}X+\delta_{2}\bar{\Delta}\bar{\Delta}X+\delta_{2}\sum_{i=1}^{m}\big{[}(\nabla_{e_{i}}R)(e_{i},S(X))X$
$\displaystyle+R(e_{i},\nabla_{e_{i}}S(X))X+2R(e_{i},S(X))\nabla_{e_{i}}X\big{]})\Big{\\}}v_{g}.$

Then, the desired formula follows. ∎

######  Remark 3.5.

As usual, Theorem 3.4 holds when $(M,g)$ is a non-compact Riemannian manifold.
Indeed, if $M$ is non-compact, we can take an open subset $W$ in $M$ whose
closure is compact, and take an arbitrary $V$ but, the support of $V$ is
contained in $W$ (see Proposition 3 of [10]).

Then, we deduce the following.

######  Corollary 3.6.

A vector field $X$ of an $m$-dimensional Riemannian manifold $(M,g)$ is
interpolating sesqui-harmonic if and only if

$\delta_{1}\bar{\Delta}X+\delta_{2}\bar{\Delta}\bar{\Delta}X+\delta_{2}\sum_{i=1}^{m}[(\nabla_{e_{i}}R)(e_{i},S(X))X\\\
+R(e_{i},\nabla_{e_{i}}S(X))X+2R(e_{i},S(X))\nabla_{e_{i}}X]=0,$

where $\\{e_{i}\\}_{i=1}^{m}$ is a local orthonormal frame field of $(M,g)$.

A reformulation of Theorem 3.2 is then

######  Corollary 3.7.

Let $(M,g)$ be an $m$-dimensional Riemannian manifold and
$X\in\mathfrak{X}(M)$. Then $X$ is an interpolating sesqui-harmonic map of
$(M,g)$ to $(TM,g_{S})$ if and only if $X$ is interpolating sesqui-harmonic
vector field and

$\delta_{1}S(X)+\delta_{2}\bar{\Delta}S(X)+\delta_{2}R(X,\bar{\Delta}X)S(X)-\delta_{2}\sum_{i=1}^{m}[R(X,\nabla_{e_{i}}\bar{\Delta}X)e_{i}\\\
-R(\nabla_{e_{i}}X,\bar{\Delta}X)e_{i}-R(e_{i},S(X))e_{i}-(\nabla_{S(X)}R)(\nabla_{e_{i}}X,X)e_{i}\\\
+R(X,\nabla_{e_{i}}X)\nabla_{e_{i}}S(X)-R(X,R(e_{i},S(X))X)e_{i}]=0.$

In the following Theorem, when $\delta_{1}$ and $\delta_{2}$ have the same
sign, we study the condition under of which a vector field $X$ of a Riemannian
manifold $(M,g)$ is interpolating sesqui-harmonic under the assumption that
the base manifold $(M,g)$ is compact. In particular, we have

######  Theorem 3.8.

Let $(M,g)$ be a compact oriented $m$-dimensional Riemannian manifold and
$X\in\mathfrak{X}(M)$ a vector field. Assuming that $\delta_{1}$ and
$\delta_{2}$ have the same sign. Then, $X:(M,g)\longrightarrow(TM,g_{S})$ is
interpolating sesqui-harmonic vector field if and only if $X$ is parallel.

###### Proof.

We assume that $X$ is an interpolating sesqui-harmonic vector field i.e.
critical point of the energy functional $E_{\delta_{1},\delta_{2}}$ restricted
to the space $\mathfrak{X}(M)$. We consider the smooth 1-parameter variation
$X_{t}=(1+t)X$ of $X$ ($|t|<\epsilon$). By the Lemma 2.15 in [7] we have

$g(\bar{\Delta}X,X)=\frac{1}{2}\Delta(\left\|X\right\|^{2})+\left\|\nabla
X\right\|^{2},$

where $\Delta$ is the Laplace-Betrami operator on functions. Applying the
divergence Theorem for the function $\left\|X\right\|^{2}$, we get

$\int_{M}g(\bar{\Delta}X,X)v_{g}=\int_{M}\left\|\nabla X\right\|^{2}v_{g}.$
(3.4)

Thus, using (3.4) and (3.4) and the proof of Theorem 3.4 in [17] we deduce
that

$\displaystyle 0$
$\displaystyle=\frac{d}{dt}E_{\delta_{1},\delta_{2}}(X_{t})\bigg{|}_{t=0}=2\delta_{1}\int_{M}g(\bar{\Delta}X,X)v_{g}+2\delta_{2}\int_{M}g(\bar{\Delta}X,\bar{\Delta}X)v_{g}$
$\displaystyle\quad+4\delta_{2}\int_{M}g(S(X),S(X))v_{g}$
$\displaystyle=2\delta_{1}\int_{M}\left\|\nabla
X\right\|^{2}v_{g}+2\delta_{2}\int_{M}\left\|\bar{\Delta}X\right\|^{2}v_{g}+4\delta_{2}\int_{M}\left\|S(X)\right\|^{2}v_{g}.$

Since $\delta_{1}$ and $\delta_{2}$ have the same sign and both functions
$\left\|\nabla X\right\|^{2}$, $\left\|\bar{\Delta}X\right\|^{2}$,
$\left\|S(X)\right\|^{2}$ are positive, we conclude that $\nabla X=0$ i.e. $X$
is parallel. Conversely, we assume that the vector field $X$ is parallel, by
virtue of Corollary 3.6, $X$ is an interpolating sesqui-harmonic vector field.
∎

######  Theorem 3.9.

Let $(M,g)$ be a compact oriented $m$-dimensional Riemannian manifold and
$X\in\mathfrak{X}(M)$ a vector field. Assuming that $\delta_{1}$ and
$\delta_{2}$ have the same sign, then $X:(M,g)\longrightarrow(TM,g_{S})$ is an
interpolating sesqui-harmonic map if and only if $X$ is parallel.

###### Proof.

We assume that $X:(M,g)\longrightarrow(TM,g_{S})$ is an interpolating sesqui-
harmonic map, then from Corollary 3.6, $X$ is an interpolating sesqui-harmonic
vector field, hence by Theorem 3.8 $X$ is parallel. Conversely, we suppose
that the vector field $X$ is parallel, by virtue of Theorem 3.2, $X$ is an
interpolating sesqui-harmonic map. ∎

######  Example 3.10.

Consider the Sol space as the Cartesian $3$-space $\mathbb{R}^{3}(x,y,z)$
equipped with the metric $g_{Sol}=e^{2z}(dx)^{2}+e^{-2z}(dy)^{2}+(dz)^{2}$ and
the orthonormal frame field $\\{e_{1}=e^{-z}\frac{\partial}{\partial
x},e_{2}=e^{z}\frac{\partial}{\partial y},e_{3}=\frac{\partial}{\partial
z}\\}.$ We consider the vector field $X=f(z)e_{3}$, where $f(z)$ is a smooth
real function depending of the variable $z$. Drawing on computations from [1],
we obtain that $X=f(z)e_{3}$ is an interpolating sesqui-harmonic map if and
only if the function $f$ satisfies the following homogeneous fourth order
differential equation.

$\delta_{2}f^{\prime\prime\prime\prime}-(\delta_{1}+4\delta_{2})f^{\prime\prime}+(2\delta_{1}+4\delta_{2})f=0.$
(3.5)

If $\frac{\delta_{1}+2\delta_{2}}{\delta_{2}}>0$, the general solution of
(3.5) is

$f(z)=c_{1}e^{\sqrt{2}z}+c_{2}e^{-\sqrt{2}z}+c_{3}e^{\sqrt{\frac{\delta_{1}+2\delta_{2}}{\delta_{2}}}z}+c_{4}e^{-\sqrt{\frac{\delta_{1}+2\delta_{2}}{\delta_{2}}}z},$
(3.6)

where $c_{1},c_{2},c_{3}$ and $c_{4}$ are real constants. In particular,
$X=f(z)e_{3}$ is also interpolating sesqui-harmonic vector field, where $f(z)$
is given by (3.6).

##  4\. Interpolating sesqui-harmonicity of vector fields on Lie groups

In this section, we investigate the interpolating sesqui-harmonicity of left-
invariant vector fields on Lie groups. We know that the action of any discrete
subgroup $\Gamma$ of a Lie group $G$ by left translations is free and properly
discontinuous. Consequently, the set of orbits, that is, the space of right
cosets $\Gamma\backslash G$, is a $C^{\infty}$-manifold and the naturel
projection $\pi:G\longrightarrow\Gamma\backslash G$, applying each $x$ to its
orbit $\Gamma x$, is a $C^{\infty}$-mapping (see [3]). Furthermore, each left-
invariant vector field on $G$ descends to $\Gamma\backslash G$, or
equivalently, if $X$ is left-invariant, then
$\pi_{\ast}X_{ba}=\pi_{\ast}X_{a}$, for all $a\in G$ and $b\in\Gamma$, we have
also each left-invariant metric on $G$ and, in general, all its left-invariant
tensor fields, descend to the quotient space $\Gamma\backslash G$ (see [11]).
Thereby, $\Gamma\backslash G$ is a Riemannian manifold with the same curvature
properties for the curvature tensor as on $G$ (see [11]). Thus, we deduce that
the projections of left-invariant vector fields preserve the properties to be
harmonic, biharmonic, interpolating sesqui-harmonic and to determine harmonic
maps, biharmonic maps and interpolating sesqui-harmonic maps.

Now, we focus on the case of compact $\Gamma\backslash G$. It should be noted
that a necessary and sufficient condition for compactness is the existence of
a compact subset $K\subset G$ whose $\Gamma$-orbits cover $G$, that is,
$\Gamma K=G$ (see [11]). In particular, any three-dimensional unimodular Lie
group $G$ admits a discrete subgroup $\Gamma$ such that $\Gamma\backslash G$
is compact and conversely (see [18]).

######  Theorem 4.1.

Let $G$ be a Lie group having a discrete subgroup $\Gamma$ such that
$\Gamma\backslash G$ is compact. Assuming that $\delta_{1}$ and $\delta_{2}$
have the same sign, then a left-invariant vector field $X$ on $G$ is
interpolating sesqui-harmonic (resp. interpolating sesqui-harmonic map) if and
only if $X$ is parallel.

###### Proof.

Suppose that $X$ is a left-invariant vector field on $G$ which is non-parallel
interpolating sesqui-harmonic, it follows that its projection on
$\Gamma\backslash G$ which is also denoted by $X$ is non-parallel
interpolating sesqui-harmonic and as $\Gamma\backslash G$ is compact, this
gives a contradiction by combining with Theorem 3.8, so we must have $X$ is
parallel. Using the similar way, we can prove the result when $X$ is an
interpolating sesqui-harmonic map by virtue of Theorem 3.9. ∎

######  Corollary 4.2.

Let $G$ be a Lie group having a discrete subgroup $\Gamma$ such that
$\Gamma\backslash G$ is compact. Then

  1. (1)

A left-invariant vector field $X$ on $G$ is harmonic (resp. harmonic map) if
and only if $X$ is parallel.

  2. (2)

A left-invariant vector field $X$ on $G$ is biharmonic (resp. biharmonic map)
if and only if $X$ is parallel.

######  Corollary 4.3.

Let $G$ be a three-dimensional unimodular Lie group. Assuming that
$\delta_{1}$ and $\delta_{2}$ have the same sign, then a left-invariant vector
field $X$ on $G$ is interpolating sesqui-harmonic (resp. interpolating sesqui-
harmonic map) if and only if $X$ is parallel.

######  Corollary 4.4.

Let $G$ be a three-dimensional unimodular Lie group. Then

  1. (1)

A left-invariant vector field $X$ on $G$ is harmonic (resp. harmonic map) if
and only if $X$ is parallel.

  2. (2)

A left-invariant vector field $X$ on $G$ is biharmonic (resp. biharmonic map)
if and only if $X$ is parallel.

######  Remark 4.5.

We have proved the Corollary 4.4 see Theorems 4.1 and 4.2 in [1].

Now, we completely determine the set of left-invariant interpolating sesqui-
harmonic vector fields and left-invariant vector fields which are
interpolating sesqui-harmonic maps on the Heisenberg group $Nil$ considered as
the Cartesian $3$-space $\mathbb{R}^{3}(x,y,z)$ equipped with the left-
invariant metric $g_{Nil}$ given by $g_{Nil}=(dx)^{2}+(dy-xdz)^{2}+(dz)^{2}.$
Since the Heisenberg group $Nil$ is three-dimensional unimodular Lie group,
then by the Corollary 4.3 we will restrict ourselves to the case where
$\delta_{1}$ and $\delta_{2}$ have different signs. The left-invariant vector
fields $\\{e_{1}=\frac{\partial}{\partial x},\;e_{2}=\frac{\partial}{\partial
y},\;e_{3}=\frac{\partial}{\partial z}+x\frac{\partial}{\partial y}\\}$
constitute an orthonormal basis of the Lie algebra $\mathfrak{g}$ of $Nil$.
The corresponding components of the Levi-Civita connection are determined by

$\displaystyle\nabla_{e_{1}}e_{1}$ $\displaystyle=0,$
$\displaystyle\nabla_{e_{1}}e_{2}$ $\displaystyle=-\frac{1}{2}e_{3},$
$\displaystyle\nabla_{e_{1}}e_{3}$ $\displaystyle=\frac{1}{2}e_{2},$
$\displaystyle\nabla_{e_{2}}e_{1}$ $\displaystyle=-\frac{1}{2}e_{3},$
$\displaystyle\nabla_{e_{2}}e_{2}$ $\displaystyle=0,$
$\displaystyle\nabla_{e_{2}}e_{3}$ $\displaystyle=\frac{1}{2}e_{1},$ (4.1)
$\displaystyle\nabla_{e_{3}}e_{1}$ $\displaystyle=-\frac{1}{2}e_{2},$
$\displaystyle\nabla_{e_{3}}e_{2}$ $\displaystyle=\frac{1}{2}e_{1},$
$\displaystyle\nabla_{e_{3}}e_{3}$ $\displaystyle=0.$

Also the curvature components are given by

$\displaystyle R(e_{1},e_{2})e_{1}$ $\displaystyle=-\frac{1}{4}e_{2},$
$\displaystyle R(e_{1},e_{2})e_{2}$ $\displaystyle=\frac{1}{4}e_{1},$
$\displaystyle R(e_{1},e_{2})e_{3}$ $\displaystyle=0,$ $\displaystyle
R(e_{2},e_{3})e_{1}$ $\displaystyle=0,$ $\displaystyle R(e_{2},e_{3})e_{2}$
$\displaystyle=-\frac{1}{4}e_{3},$ $\displaystyle R(e_{2},e_{3})e_{3}$
$\displaystyle=\frac{1}{4}e_{2},$ (4.2) $\displaystyle R(e_{3},e_{1})e_{1}$
$\displaystyle=-\frac{3}{4}e_{3},$ $\displaystyle R(e_{3},e_{1})e_{2}$
$\displaystyle=0,$ $\displaystyle R(e_{3},e_{1})e_{3}$
$\displaystyle=\frac{3}{4}e_{1}.$

Let $X=\alpha e_{1}+\beta e_{2}+\gamma e_{3}$ an arbitrary left-invariant
vector field on $(Nil,g_{Nil})$. By using (4) and (4), we yield

$\displaystyle\bar{\Delta}X$
$\displaystyle=\frac{\alpha}{2}e_{1}+\frac{\beta}{2}e_{2}+\frac{\gamma}{2}e_{3},$
$\displaystyle\bar{\Delta}\bar{\Delta}X$
$\displaystyle=\frac{\alpha}{4}e_{1}+\frac{\beta}{4}e_{2}+\frac{\gamma}{4}e_{3},$
$\displaystyle S(X)$
$\displaystyle=\frac{-\beta\gamma}{4}e_{1}+\frac{\alpha\beta}{4}e_{3}.$ (4.3)

By virtue of (4)-(4), a long but straightforward calculation we get

$\displaystyle\delta_{1}$
$\displaystyle\bar{\Delta}X+\delta_{2}\bar{\Delta}\bar{\Delta}X+\delta_{2}\sum_{i=1}^{3}[(\nabla_{e_{i}}R)(e_{i},S(X))X+R(e_{i},\nabla_{e_{i}}S(X))X$
$\displaystyle+2R(e_{i},S(X))\nabla_{e_{i}}X]=\delta_{1}\bar{\Delta}X+\delta_{2}\bar{\Delta}\bar{\Delta}X+\delta_{2}\sum_{i=1}^{3}[\nabla_{e_{i}}R(e_{i},S(X))X$
$\displaystyle+R(e_{i},S(X))\nabla_{e_{i}}X]=\frac{\alpha(8\delta_{1}+\delta_{2}(4+\beta^{2}))}{16}e_{1}$
$\displaystyle+\frac{\beta(8\delta_{1}+\delta_{2}(4+\alpha^{2}+\gamma^{2}))}{16}e_{2}+\frac{\gamma(8\delta_{1}+\delta_{2}(4+\beta^{2}))}{16}e_{3}.$
(4.4)

On the other hand using (4)-(4), a long but direct calculations we find

$\bar{\Delta}S(X)=-\frac{\beta\gamma}{8}e_{1}+\frac{\alpha\beta}{8}e_{3},$
(4.5) $R(X,\bar{\Delta}X)S(X)=0,$ (4.6)
$\sum_{i=1}^{3}R(X,\nabla_{e_{i}}\bar{\Delta}X)e_{i}=\frac{\beta\gamma}{8}e_{1}-\frac{\alpha\beta}{8}e_{3},$
(4.7)
$\sum_{i=1}^{3}R(\nabla_{e_{i}}X,\bar{\Delta}X)e_{i}=-\frac{\beta\gamma}{8}e_{1}+\frac{\alpha\beta}{8}e_{3},$
(4.8)
$\sum_{i=1}^{3}R(e_{i},S(X))e_{i}=-\frac{\beta\gamma}{8}e_{1}+\frac{\alpha\beta}{8}e_{3},$
(4.9)
$\sum_{i=1}^{3}(\nabla_{S(X)}R)(\nabla_{e_{i}}X,X)e_{i}=\frac{\beta\gamma(\alpha^{2}+\beta^{2}+\gamma^{2})}{16}e_{1}-\frac{\alpha\beta(\alpha^{2}+\beta^{2}+\gamma^{2})}{16}e_{3},$
(4.10)
$\sum_{i=1}^{3}R(X,\nabla_{e_{i}}X)\nabla_{e_{i}}S(X)=\frac{\beta\gamma(3\alpha^{2}+3\gamma^{2}-\beta^{2})}{64}e_{1}-\frac{\alpha\beta(3\alpha^{2}+3\gamma^{2}-\beta^{2})}{64}e_{3},$
(4.11)
$\sum_{i=1}^{3}R(X,R(e_{i},S(X))X)e_{i}=-\frac{\beta\gamma(9\alpha^{2}+9\gamma^{2}+\beta^{2})}{64}e_{1}+\frac{\alpha\beta(9\alpha^{2}+9\gamma^{2}+\beta^{2})}{64}e_{3}.$
(4.12)

From (4.5)-(4.12), we yield

$\displaystyle\delta_{1}S(X)+\delta_{2}\bar{\Delta}S(X)+\delta_{2}R(X,\bar{\Delta}X)S(X)-\delta_{2}\sum_{i=1}^{3}[R(X,\nabla_{e_{i}}\bar{\Delta}X)e_{i}$
$\displaystyle-R(\nabla_{e_{i}}X,\bar{\Delta}X)e_{i}-R(e_{i},S(X))e_{i}-(\nabla_{S(X)}R)(\nabla_{e_{i}}X,X)e_{i}$
$\displaystyle+R(X,\nabla_{e_{i}}X)\nabla_{e_{i}}S(X)-R(X,R(e_{i},S(X))X)e_{i}]$
$\displaystyle=-\frac{\beta\gamma(4\delta_{1}+\delta_{2}(8+\alpha^{2}+\gamma^{2}-2\beta^{2}))}{16}e_{1}+\frac{\alpha\beta(4\delta_{1}+\delta_{2}(8+\alpha^{2}+\gamma^{2}-2\beta^{2}))}{16}e_{3}.$
(4.13)

From (4) and (4), we deduce that $X$ is an interpolating sesqui-harmonic map
if and only if

$\left\\{\begin{array}[]{ll}\alpha(8\delta_{1}+\delta_{2}(4+\beta^{2}))=0,&\hbox{}\\\
\beta(8\delta_{1}+\delta_{2}(4+\alpha^{2}+\gamma^{2}))=0,&\hbox{}\\\
\gamma(8\delta_{1}+\delta_{2}(4+\beta^{2}))=0,&\hbox{}\end{array}\right.$
(4.14)

and

$\left\\{\begin{array}[]{ll}\beta\gamma(4\delta_{1}+\delta_{2}(8+\alpha^{2}+\gamma^{2}-2\beta^{2}))=0,&\hbox{}\\\
\alpha\beta(4\delta_{1}+\delta_{2}(8+\alpha^{2}+\gamma^{2}-2\beta^{2}))=0.&\hbox{}\end{array}\right.$
(4.15)

In particular, $X$ is an interpolating sesqui-harmonic vector field if and
only if (4.14) holds. From the system (4.14), the subcases $\beta=\gamma=0$,
$\alpha=\gamma=0$ and $\alpha=\beta=0$ give vector fields $X=\alpha e_{1}$,
$X=\beta e_{2}$ and $X=\gamma e_{3}$ respectively such that
$\delta_{1}=-\frac{1}{2}\delta_{2}$ which define by (4.15) an interpolating
sesqui-harmonic maps. If $\alpha=0$ (resp. $\gamma=0$), the system (4.14) give
the solutions $X=\pm
2\sqrt{-\frac{(2\delta_{1}+\delta_{2})}{\delta_{2}}}e_{2}\pm
2\sqrt{-\frac{(2\delta_{1}+\delta_{2})}{\delta_{2}}}e_{3}$ (resp. $X=\pm
2\sqrt{-\frac{(2\delta_{1}+\delta_{2})}{\delta_{2}}}e_{1}\pm
2\sqrt{-\frac{(2\delta_{1}+\delta_{2})}{\delta_{2}}}e_{2}$) whenever
$\frac{(2\delta_{1}+\delta_{2})}{\delta_{2}}<0$, which define by (4.15) an
interpolating sesqui-harmonic maps if $\delta_{1}=-\delta_{2}$. If $\beta=0$,
the system (4.14) gives the solution $X=\alpha e_{1}+\gamma e_{3}$ such that
$\delta_{1}=-\frac{1}{2}\delta_{2}$ which defines by (4.15) an interpolating
sesqui-harmonic map. If $\alpha\neq 0$, $\beta\neq 0$, $\gamma\neq 0$, the
system (4.14) admits $\beta=\pm
2\sqrt{-\frac{(2\delta_{1}+\delta_{2})}{\delta_{2}}}$ and the cylinder
$\alpha^{2}+\gamma^{2}=-4\frac{(2\delta_{1}+\delta_{2})}{\delta_{2}}$ whenever
$\frac{(2\delta_{1}+\delta_{2})}{\delta_{2}}<0$. Therefore, the coordinates of
$X$ satisfy the equations of circles:

$\displaystyle C^{1}$
$\displaystyle=\biggl{\\{}\alpha^{2}+\gamma^{2}=-4\frac{(2\delta_{1}+\delta_{2})}{\delta_{2}},\;\beta=2\sqrt{-\frac{(2\delta_{1}+\delta_{2})}{\delta_{2}}}\biggr{\\}};$
$\displaystyle C^{2}$
$\displaystyle=\biggl{\\{}\alpha^{2}+\gamma^{2}=-4\frac{(2\delta_{1}+\delta_{2})}{\delta_{2}},\;\beta=-2\sqrt{-\frac{(2\delta_{1}+\delta_{2})}{\delta_{2}}}\biggr{\\}},$

and by virtue of (4.15) the corresponding vector fields define an
interpolating sesqui-harmonic maps if $\delta_{1}=-\delta_{2}$.

Summarizing, we yield

######  Theorem 4.6.

Let $X=\alpha e_{1}+\beta e_{2}+\gamma e_{3}$ be a left-invariant vector field
on the Heisenberg group $(Nil,g_{Nil})$. Then $X$ is an interpolating sesqui-
harmonic vector field if and only if

  * •

$X=\alpha e_{1},\;\beta e_{2},\;\gamma e_{3},$ with
$\delta_{1}=-\frac{1}{2}\delta_{2}$. The corresponding vector fields define an
interpolating sesqui-harmonic maps.

  * •

$X=\pm 2\sqrt{-\frac{(2\delta_{1}+\delta_{2})}{\delta_{2}}}e_{2}\pm
2\sqrt{-\frac{(2\delta_{1}+\delta_{2})}{\delta_{2}}}e_{3}$ (resp. $X=\pm
2\sqrt{-\frac{(2\delta_{1}+\delta_{2})}{\delta_{2}}}e_{1}\pm
2\sqrt{-\frac{(2\delta_{1}+\delta_{2})}{\delta_{2}}}e_{2}$) whenever
$\frac{(2\delta_{1}+\delta_{2})}{\delta_{2}}<0$, which define an interpolating
sesqui-harmonic maps if $\delta_{1}=-\delta_{2}$.

  * •

$X=\alpha e_{1}+\beta e_{2}+\gamma e_{3}$ such that the coordinates of $X$
satisfy the equations of circles $C^{1}$ and $C^{2}$. The corresponding vector
fields define an interpolating sesqui-harmonic maps if
$\delta_{1}=-\delta_{2}$.

######  Remark 4.7.

It should be noted already that a vector field is interpolating sesqui-
harmonic does not automatically imply that the corresponding map is
interpolating sesqui-harmonic map. Indeed, from Theorem 4.6 if
$\delta_{1}\neq-\delta_{2}$, the vector fields $X=\pm
2\sqrt{-\frac{(2\delta_{1}+\delta_{2})}{\delta_{2}}}e_{1}\pm
2\sqrt{-\frac{(2\delta_{1}+\delta_{2})}{\delta_{2}}}e_{2}$ are interpolating
sesqui-harmonic but do not define an interpolating sesqui-harmonic maps.

## References

  * [1] Alem, A., Kacimi, B., Özkan, M.: Vector fields which are biharmonic maps. J. of geometry. 113 14, (2022) DOI https://doi.org/10.1007/s00022-022-00627-5.
  * [2] Baird, P., Wood, J.C.: Harmonic morphisms between Riemannian manifolds. In: Lond. Math. Soc. Monogr., vol. 29. Oxford University Press, Oxford (2003)
  * [3] Boothby, W.M.: An Introduction to Differentiable Manifolds and Riemannian Geometry. Academic Press, New York (1986)
  * [4] Branding, V.: On interpolating sesqui-harmonic maps between riemannian manifolds. J. Geom. Anal. 30, 248–273 (2020).
  * [5] Branding, V.: Some analytic results on interpolating sesqui-harmonic maps. Annali. di. Matematica. 199, 2039–2059 (2020).
  * [6] Dombrowski, P: On the geometry of tangent bundle. J. Reine Angew. Math. 210, 73–88 (1962).
  * [7] Dragomir, S., Perrone, D.: Harmonic Vector fields: Variational Principles and Differential Geometry. Elsevier, Amsterdam (2011)
  * [8] Eells, J., Lemaire, L.: Selected topics in harmonic maps, in: CBMS Regional Conference Series in Mathematics, vol. 50, Amer. Math. Soc., Providence, 1983.
  * [9] Eells, J., Sampson, J.H.: Harmonic mappings of Riemannian manifolds. Amer. J. Math. 86(1), 109–160 (1964).
  * [10] Gil-Medrano, O.: Relationship between volume and energy of unit vector fields. Differ. Geo. Appl. 15, 137–152 (2001).
  * [11] González-Dávila, J.C., Vanhecke, L.: Energy and volume of unit vector fields on three-dimensional Riemannian manifolds. Differ. Geo. Appl. 16, 225–244 (2002).
  * [12] Hélein, F., Wood, J.C.: Harmonic maps. In: Handbook of global analysis, pp 417–491, 1213. Elsevier Sci. B. V., Amsterdam (2008)
  * [13] Ishihara, T.: Harmonic sections of tangent bundles. J. Math. Univ. Tokushima. 13, 23–27 (1979).
  * [14] Jiang, G.: 2-Harmonic maps and their first and second variational formulas. Translated into English by Hajime Urakawa. Note Mat. 28(suppl. n. 1), 209–232 (2008).
  * [15] Karaca, F., Özgür, C., De, U.C: On interpolating sesqui-harmonic Legendre curves in Sasakian space forms. Int J. Geom. Methods Mod. Phys. 17(1), 2050005 (2020).
  * [16] Kowalski, O.: Curvature of the induced Riemannian metric on the tangent bundle of a Riemannian manifold. J. Reine Angew. Math. 250, 124–129 (1971).
  * [17] Markellos, M., Urakawa, H.: The biharmonicity of sections of the tangent bundle. Monatsh. Math. 178, 389–404 (2015).
  * [18] Milnor, J.: Curvature of left invariant metrics on Lie groups. Adv. Math. 21, 293–329 (1976).
  * [19] Nouhaud, O.: Applications harmoniques d’une variété Riemannienne dans son fibré tangent. C R Acad. Sci. Paris 284, 815–818 (1977).
  * [20] Ou, Y.L, Chen, B.Y.: Biharmonic Submanifolds and Biharmonic Maps in Riemannian Geometry. World Scientific (2020)
  * [21] Polyakov, A.: Fine structure of strings, Nucl. Phys. B 268(2), 406–412 (1986).
  * [22] Urakawa, H.: Geometry of Biharmonic Mappings: Differential Geometry of Variational Methods. World Scientific (2018)
  * [23] Yano, K., Ishihara, S.: Tangent and Cotangent Bundles Dekker, New York (1973)