# Investigation of correlation effects in FeSe superconductor by LDA+U method

H. Lohani, P. Mishra and B. R. Sekhar

###### Abstract

Correlation effects are observed strong in Iron chalcogenides superconductors
by experimental and theoretical investigations. We present a comparative study
of the influence of Coulomb interaction and Hund’s coupling in the electronic
structure of FeSe and FeTe. The calculation is based on density functional
theory (DFT) with local density approximation(LDA+U) framework employed in TB-
LMTO ASA code. We found the correlation effects were orbital selective due to
the strength of interorbital hybridization among different Fe-3d orbitals
mediated via chalcogen (Se/Te-p) orbitals is different in both the compounds,
however Coulomb interaction is screened significantly by Te-p bands in FeTe.
Similarly the orbital section is different in both the compounds because of
the difference in the chalcogen height.

##  1 Introduction

Iron based superconductors, particularly members of FeSe1-xTex family attract
much attention due to their nature of strong electron correlation unlike other
superconductors. A recent advancement in this field is the synthesis of
single-layer films of FeSe on SrTiO3 substrates exhibiting superconductivity
(Tc=80 K) which turns insulating with the addition of one more layer[1]. This
unusual behavioral difference between single and double layer films of FeSe is
a signature of strong electron correlation which has been experimentally
observed[2, 3]. Superconductivity in the FeSe1-xTex compounds was first
reported by Hsu et al. [4] in the FeSe (x = 0) compound exhibiting a Tc around
8 K which rises up to 37 K under pressure (7GPa)[5]. On the other hand, the
other extreme composition of this family, Fe1.068Te, though not a
superconductor show a spin density wave (SDW) ordering at 67K[12] with an
accompanying structural transition from tetragonal to monoclinic. With Se
doping superconductivity emerges in FeTe with a simultaneous decrease in the
SDW and the value of Tc reaches maximum 15 K for x = 0.5 doped case. The Fe
content is also detrimental for superconductivity; with excess Fe favors the
spin localization destroying the superconductivity[14, 15]. Both FeSe and FeTe
have tetragonal crystal structure belonging to space group symmetry P4/nmm. It
consists a square planar sheet of Fe atoms, which is tetrahedraly coordinated
with anion (Se/Te) atoms. However, the height of anion atom from the Fe square
plane is different in these two compounds and this plays a pivotal role in
determining the electronic properties of these systems[10, 13].

A recent ARPES study on FeSe1-xTex compositions by Ieki et al. [39] has shown
clearly the strong electronic correlation in these compounds, where a small
quasi particle weight in FeTe transforms into a sharp one with increase in Se
content. Other ARPES results [32, 37] on these compounds have shown
significant band renormalization which was supported by the calculations based
on local density approximation (LDA). Tamai et al. [35] has found the mass
renormalization factor to be m∗/m = 20 from their ARPES study on
FeSe0.42Te0.58. It is close to the value observed in highly correlated systems
like transition metal oxides. Also, our angle integrated valence band
photoemission study on FeSe1-xTex [19] revealed significant spectral weight
shifts in the near Ef region with Se doping leading to the formation of a
pseudogap. Further, a temperature dependent orbital selective spectral weight
transfer was also reported by us[19]. Although, such manifestations of the
strong coulomb correlation were also shown in many other photoemisson
studies[20, 36], these experimental observations are not well addressed by LDA
based electronic structure calculations. However, results of some recent
calculations[17, 18, 31], where Coulomb correlations were included by using
LDA+DMFT frame work, are very close to the experimental findings. In this
paper we are presenting our calculations showing the evolution of electronic
structure by the incorporation of different values of Coulomb interaction U
and intra-atomic exchange J based on LDA+U scheme in FeSe and FeTe. We
observed multi orbital correlation effect in Fe-3d states which is more
prominent in FeSe in comparison to FeTe. We have discussed our results
referring to the difference in the geometry of the anion tetrahedra in both
the compounds.

##  2 Details of Calculation

The band structure calculations we performed were based on first principles
Langreth-Meh-Hu gradient corrected von Barth Hedin parametrized LDA [27]
energy and potential. Lattice parameters used in the calculations are taken
from the experimental data[21, 22] published earlier by others. Correlation
effects of the Fe-3d orbitals have been examined by employing different values
of Coulomb interaction parameter U and Hund’s coupling J[28]. Empty spheres
were introduced to make the volume of the unit cell equal to the total volume
of the spheres within the permissible limit of atomic sphere approximation
(ASA). Fe-4s, 4p, 3d ; Se-4s, 4p, 4d and Te-5s, 5p, 5d orbitals are used as
the basis set for the valence energy region. A mesh of 12$\times$12$\times$8
is used for sampling the irreducible part of the Brillouin zone integration.
Height of the anion atom from the Fe square plane was relaxed by minimizing
the total energy using Quantum Espresso code[29]. We checked our calculation
parameters by comparing them with the LDA results reported earlier[36, 30].

##  3 Results and Discussion

###  3.1 DOS

Fig (1a) and (1b) show the plots of density of states (DOS) for FeSe and FeTe
respectively for the range of 4.0 to -7.0 eV binding energy (BE). The valence
band (VB) DOS comprises of the region from Ef to -6.0 eV BE while the
unoccupied DOS extends from Ef to 3.0 eV BE. The near Ef states from 2.8 to
-2.4 eV are predominantly Fe-3d derived for both FeSe and FeTe. These Fe-3d
states, separated in lower and upper bands exhibit a clear pseudogap feature
just above the Ef (0.24 eV) in the case of FeSe whereas it is less prominent
in the case of FeTe. The states around -2.2 to -6.0 eV are originate from the
Fe-3d and anion (Se/Te)-p hybridized states. Interestingly, for the case of
FeSe, there exists a sharp gap at -2.3 eV which is not present in case of
FeTe. Due to the smaller electronegativity of Te compared to Se, the
hybridized states between Fe-3d and anion-p orbitals are placed at a lower BE
in FeTe in comparision to those of FeSe. DOS have also been calculated after
downfolding the valence orbital of the anion atom (Fig (1c) and (1d)) in order
to get a deeper insight into the role of anion orbitals. The gap at -2.3 eV
present in Fig (1a) is due to the splitting of bonding and antibonding bands
between Fe-3d and Se-p hybridized states which becomes less prominent when the
valence orbitals of Se atom are downfolded as shown in Fig (1d). Similarly the
pseudogap feature across the Ef is also associated with the hybridization
between Fe-3d and Se-p orbitals. Unlike FeSe, the anion orbitals do not play
any major role in modifying the DOS in FeTe as is clear from Fig (1b) and
(1d). It indicates a weak hybridization between Fe-3d and Te-p orbitals in
FeTe. The role of the anion orbitals is linked to the structural geometry of
FeSe and FeTe. The insets of Fig (1c) and (1d) show the geometry of anion
tetrahedra in FeSe and FeTe respectively. This tetrahedral geometry depends on
two important parameters; firstly, the height of anion from the Fe square
plane (z) which is 1.64 Å and 1.46 Å in case of FeTe and FeSe respectively and
secondly, the anion-Fe-anion angle ($\alpha$). The enhanced z height in case
of FeTe reduces the interorbital hoping among the Fe-3d orbitals, mediated via
anion p orbitals. Similarly, the value of $\alpha$ is 99.9∘ in case of FeTe
which increases to the perfect tetrahedron value 109.4∘ in case of FeSe. The
large value of $\alpha$ and small value of anion height (z) makes a stronger
hybridization between the Fe-3d and anion-p orbitals in case of FeSe in
comparison to FeTe. This difference in hybridization strength is reflected in
the plot of DOS Fig (1a) and Fig (1c).

Coulomb correlation effects are important in bands of narrow width, especially
in Fe-3d states. So the changes in Fe-3d states under the influence of
different values of U have been calculated and shown for FeSe and FeTe
respectively in Panel (a) and (b) of Fig 2. In FeSe, Fe-3d states start
localizing with the application of U and noticeable changes occur at higher
values of U. For U = 4.0 eV case, the pseudogap feature disappears around Ef
and only two peaks are observed in the VB comparared to the case of U = 0.0
eV. These two peaks merge and shift towards higher BE with U = 5.0 eV. In case
of FeTe, the Fe-3d states also become localized under the application of U.
However the amount of shift towards higher BE at large values of U (3.5 and
5.0 eV) is less and states are more at the vicinity of Ef at smaller values of
U (1.0 and 2.0 eV). In addition to this, narrowing of the lower and upper
bands of Fe-3d states with increase of U is less. It indicates that the effect
of Coulomb correlation is weak in FeTe. Possible reason for this, is the
presence of Te-p states at lower BE which screen the U strongly [23]. On the
basis of previous reports [23, 17, 31] the value of U=4.0 and 3.5 eV are
chosen to see the evolution of Fe-3d states under the influence of J in FeSe
and FeTe as shown in Fig 3(a) and 3(b) respectively. It is observed that the
Fe-3d states are modified significantly even by introduction of a small value
of J=0.1 eV in case of FeSe. The Hund’s coupling, shifts all the Fe-3d states
towards lower BE, with a simultaneous appearance of pseudogap slightly above
Ef. With increase in the value of J further, there is no substantial changes
in the DOS. In the case of FeTe, the Fe-3d states are also shifted towards
lower BE, particularly the states near Ef increase gradually, with the
incorporation of J though the changes are less as compared to FeSe. These
results show that Hund’s coupling J is a key factor in the formation of Fe-3d
DOS.

In order to heighlight the correlation effects, DOS of different Fe-3d
orbitals are plotted at different values of U and J in FeSe and FeTe as shown
in Panel (a) and (b) of Fig 4 respectively. In FeSe, in the absence of U and
J, near Ef states and pseudogap feature arises from dyz/xz and
d${}_{x^{2}-y^{2}}$ orbitals. The states originating from dxy orbital have
largest splitting with two peaks at -1.7 and 1.6 eV BE in the DOS.
Additionally, a clear gap is observed in d${}_{3z^{2}-r^{2}}$ states around Ef
which are quite localized at -0.8 eV. Application of Coulomb interaction(U=4.0
eV), results in localization of the states derived from all four orbitals and
the states shift towards higher BE. This shift is the largest in the states of
d${}_{x^{2}-y^{2}}$ orbital. Major effect of U is observed in dyz/xz states,
where broad states near Ef transforms into two clear peaks at higher BE. Hence
the pseudogap feature vanishes across Ef. Application of a small value of
Hund’s coupling J=0.1 eV, restore the dyz/xz states near Ef and no significant
changes are seen by further increasing the J value from 0.1 to 1.2 eV. This
nature of pseudogap, which occurs in full range of Hund’s coupling but absent
when J=0.0 at U=4.0 eV, is consistent with the previous work of Ansgar et al.
[18] where this pseudogap is attributed to a resonance in self energy caused
by spin fluctuations. In case of FeTe in the absence of U and J values, near
Ef states and pseudogap are also formed by dyz/xz and d${}_{x^{2}-y^{2}}$
states like FeSe but a gap is present at -0.7 eV in dyz/xz states as well as
number of these states are more across the Ef. The d${}_{3z^{2}-r^{2}}$ and
d${}_{x^{2}-y^{2}}$ states shift towards higher BE by switching on the U=3.5
eV and the amount of shift is quite small in comparison to FeSe. These states
shows an incremental shift towards Ef with increase the J value in this case.
These changes in DOS are presented in table 1, where the occupancy of
electrons in different orbitals of Fe-d is tabulated, for quantitative
analysis. The occupancy of dyz/xz and d${}_{x^{2}-y^{2}}$ orbitals show an
opposite behaviour for FeSe and FeTe under the influence of Coulomb
correlation energy. On the other hand, a remarkable enhancement is observed in
the occupancy of d${}_{x^{2}-y^{2}}$ orbital after introducing j= 0.1 eV in
FeSe. It is almost twice in comparison to J=0.0 eV case. The change in the
occupancy of Fe-3d orbitals due to the effect of U and J presented here
summaries the orbital selective effect in Fe-3d orbitals in FeSe and FeTe. The
individual occupancy of four orbitals are different which are further enhanced
by Coulomb interaction and Hund’s coupling. It turns out to a orbital
selective nature of the correlation effect which is crucially depends on
individual band filling factor[43].

###  3.2 Band

Fig 5 shows the near Ef band structure of FeSe with U = 0.0 (5a), U = 4.0 eV
(5b) and FeTe with U = 0.0 (5c) and U = 3.5 eV (5d). Intially when U and J =
0.0 eV, three hole like and two electron like bands are observed at $\Gamma$
and M point respectively in both FeSe and FeTe. However, the outer hole like
bands are quasidegenerate in case of FeSe. In order to show the contribution
of different Fe-3d orbitals, the fat bands are calculated for both FeSe and
FeTe. These bands are plotted in Fig 5(e-h) for the FeSe case with U and J =
0.0 eV. Fatness of bands indicate that the innermost hole like band has a
dyz/xz and the outer two have dxy and d${}_{x^{2}-y^{2}}$ orbital characters
in both the compounds. Similarly the electron like bands are composed mainly
of d${}_{x^{2}-y^{2}}$ with little contribution from the dyz/xz orbital in
FeSe while inner electron like band (-0.23 eV) has dyz/xz and outer one (-0.47
eV) has d${}_{x^{2}-y^{2}}$ orbital character in FeTe. Another major
difference is that the Te-p bands are intermixed with Fe-d bands opposite to
that of FeSe, where Fe-d bands are quite separated from the Se-p bands. These
Te-p bands screen the effect of Coulomb interaction U in FeTe, hence the value
of U is smaller in FeTe in comparison to FeSe. After applying Coulomb
interaction degeneracy of the dxy and d${}_{x^{2}-y^{2}}$ hole like bands at
the $\Gamma$ point has been lifted as shown in Fig(5b). The
d${}_{x^{2}-y^{2}}$ band moves towards lower BE and other two dxy and dyz/xz
band move towards higher BE at $\Gamma$ point. The same band moves in opposite
direction in case of FeTe by the application of U= 3.5 eV. On the other hand,
the separation between the two electron like bands at the M point is enhanced
in FeTe unlike the FeSe case where both the electron like bands are
quasidegenerate and shift towards the lower BE under the influence of U = 4.0
eV.

Evolution of the near Fermi bands in FeSe and FeTe, with the tuning of
different values of J, are shown in panel a and b of Fig 6 respectively. In
FeSe, at $\Gamma$ point, The d${}_{x^{2}-y^{2}}$ band which was crossing the
Ef comes down below the Ef and the dyz/xz band which was not crossing Ef,
crosses Ef after applying J = 0.1 eV. Similarly, degeneracy of the two
electron like bands at M point is also lifted. In FeTe, a gradual shifting is
observed in the hole like d${}_{x^{2}-y^{2}}$ band at the $\Gamma$ point
towards the Ef and a gradual decrease in the separation between two electron
like bands at the M point under the influence of Hund’s coupling.

In a recent ARPES report on FeSe two hole like bands have been observed around
the $\Gamma$ point and two bands, one electron like and the other hole like,
at the M point from 40 meV below Ef [16]. Another ARPES study on FeTe1-xSex
for x=0, 0.2, 0.3, 0.4 and 0.45 compounds by Ieki et al. [39] reported three
clear hole like bands at the $\Gamma$ point ($\alpha$, $\alpha^{\prime}$ and
$\beta$) which evolve with Se doping and a shallow electron pocket at the M
point for x = 0.45. The $\alpha^{\prime}$ and $\beta$ crosses Ef whereas
$\alpha$ is around 20 meV below Ef. Similar results have also been
experimentally observed in Fe1.04Te0.66Se0.34 [37], FeSe0.42Te0.58 [35],
Fe1.03Te0.7Se0.3 [42] and FeTe0.55Se0.45 [40]. On the other hand in Fe1.02Te
[48] and Fe1-xTe/Se [32] only two hole like bands are observed at $\Gamma$
point. The band renormalization factor also varies highly at different points
of Brillouin zone. For example, in case of FeSe0.42Te0.58 Tamai et al. [35]
observed m*/m=20 for electron like band at the M point whereas it is just 6
for one of the hole like band at the $\Gamma$ point. It is clear from Fig (5a)
and (5c) that initially when U ans J are not incorporated all the three hole
like bands crosses the Ef in both the compounds as well as electron like band
is placed at higher BE (-0.29 eV) in FeSe and (-0.23 and -0.47 eV) in FeTe at
the M point. Only When correlation effect taken into account, one of the hole
like band comes down below the Ef and other two crosses the Ef as well as
electron like band approaches towards Ef as clear from Fig 6 where in case of
FeTe at U = 3.5 and J = 0.8 eV, innermost hole like band appear at 0.18 eV
below Ef at the $\Gamma$ point and electron like band position closer to Ef by
0.1 eV in comparison to U and J = 0.0 eV case. This trend qualitatively
matches with above mentioned experimental findings, although there is a
difference in the absolute energy scale. The orbital character of these near
Ef bands, revealed by photoemission studies using polarized light source, [37,
40] also agree with our results. These results signify the importance of
electronic correlation in these compounds however further calculation is
required to handle the correlation effects in a better way to bring down the
difference between the experimental findings and calculated results based on
LDA+U scheme.

From the transfer integral values, calculated by Miyake et al. for FeSe and
FeTe [23], it is clear that this value is large between dxy and nearest
neighbour (nn) d${}_{3z^{2}-r^{2}}$ orbital. It builds a strong interorbital
hybridization between dxy and d${}_{3z^{2}-r^{2}}$ orbitals which leads to a
localization of d${}_{3z^{2}-r^{2}}$ states at higher BE with a clear gap from
Ef. The dxy orbitals point towards the nn Fe site so it has a largest band
width and shows two well separated peaks one in the valence band and the other
in the conduction band. On the other hand, transfer integrals between dxy and
nn d${}_{x^{2}-y^{2}}$ orbital is small but for second nn it depends on the
height of anion. It is large in case of FeSe in comparison to FeTe. Thus
chalcogen-p enhance the interorbital hybridization between dxy and
d${}_{x^{2}-y^{2}}$ orbital which reflects in a prominent pseudogap structure
in d${}_{x^{2}-y^{2}}$ states in FeSe. Moreover, interorbital hybridization
between dxy and second nn dyz/xz orbital is also mediated via anion-p orbital,
hence it is also large in case of FeSe and contribute to pseudogap feature. On
the contrary, larger height of Te anion allows a finite value of transfer
integrals between dyz/xz with nn d${}_{3z^{2}-r^{2}}$ and d${}_{x^{2}-y^{2}}$
orbitals by breaking the mirror plane symmetry in case of FeTe. These inter
orbital hybridizations are responsible for the gap in dyz/xz states at -0.8
eV, which is absent in FeSe (Fig(4)). When Coulomb interaction is introduced,
it reduces the interorbital hopping and mainly dyz/xz and d${}_{x^{2}-y^{2}}$
states are affected which have large number of states near Ef. As a
consequence of this, the electron is transferred from the in-plane
d${}_{x^{2}-y^{2}}$ and dxy to the out of plane dyz/xz and
d${}_{3z^{2}-r^{2}}$ orbitals and localized them at higher BE in case of FeSe.
In FeTe, the transfer occurs from dxy and dyz/xz to d${}_{x^{2}-y^{2}}$
orbitals as is clear from Fig 2 and Fig 4. On the other hand, the application
of J blocks the fluctuations in the occupancy of different Fe-3d orbitals
which is clearly seen from the occupation table where a small value of J = 0.1
eV redistribute the electrons among different d-orbitals. A clear orbital
selective effects (increase in the occupation of dyz/xz and
d${}_{x^{2}-y^{2}}$) is seen under the influence of Hund’s coupling. Since
crystal field splitting is large in FeTe because the value of $\alpha$
deviates largely from the ideal tetrahedron value unlike the FeSe [45]. It
increases the energy difference between the d${}_{x^{2}-y^{2}}$ with the
dyz/xz and dxy orbitals. So Hund’s coupling promotes a gradual transfer of
electrons from the highly occupied d${}_{x^{2}-y^{2}}$ (at J = 0.0 eV) to
dxz/yz and dxy orbitals with an increase of J value contrary to FeSe, where a
small value of J is sufficient to transfer the electrons among these orbitals
in order to reduce their Coulomb repulsion energy. Thus the different strength
of interorbital hybridization and crystal field splitting, which is mainly
governed by the anion height, change the occupancy of electrons and band
structure of individual Fe-3d orbital. Hund’s coupling promotes this
differentiation and act like a band decoupler which was previously studied by
Medici et al. [46, 47]. This could be the origin of the orbital selective
correlation effects seen in iron chalcogenide compounds.

##  4 Conclusion

We presented a systematic study of the effect of Coulomb interaction and
Hund’s coupling in Fe-3d states in FeSe and FeTe. In both the compounds states
around Ef are predominantly originated from Fe-3d orbital having a pseudogap
feature just above the Ef, whereas hybridized states between Fe-3d and
chalcogen-p orbitals lie at higher BE. This hybridization crucially depends on
the chalcogen height from the Fe plane and it is weak in case of FeTe where
the height of Te anion is higher in comparison to Se anion height in FeSe.
Coulomb interaction localizes and shifts the Fe-3d states towards higher BE
energy in both the compounds, however this interaction is strongly screened by
Te-p bands in FeTe. It is observed that this effect is significant in dyz/xz
and d${}_{x^{2}-y^{2}}$ states in case of FeSe. Electrons in these localized
states again become itinerant under the influence of J and a clear orbital
selective changes are seen in the electronic structure. Similar to U, Hund’s
coupling effect is also prominent in FeSe in comparison to FeTe. The orbital
selective nature of the correlation effect is linked to the different values
of the interorbital hybridization among different Fe-d orbitals which is
mediated via chalcogen-p orbitals. The strength of these interorbital
hybridization mainly governs by the geometry of anion tetrahedra, height of
anion from the Fe plane (z) and anion-Fe-anion angle $\alpha$. The difference
in the anion tetrahedra geometry turns out to a different orbital selective
nature of the correlation effect in both the compounds.

## References

  * [1] Wang, Q. Y. Chin. Phys. Lett. 29 (2012) 037402.
  * [2] J. J. Lee, F. T. Schmitt, R. G. Moore, S. Johnston, Y.-T. Chi, W Li, M. Yi, Z. K. Liu, M. Hashimoto, Y. Zhang, D. H. Lu, T. P. Devereaux, D. -H. Lee and Z.-X Shen arXiv:1312.2633V5
  * [3] Xu Liu, Def Liu, W. Zhang, J. He, L. Zhao, S. He, D. Mou, F. Li et al.arXiv:1402.1400v1
  * [4] F.C.Hsu, J.-Y.Luo,K. W.Yeh, T. K. Chen, T. W. Hung, P. M. Wu, Y. C. Lee, Y. L. Huang, Y. Y. Chu, D. C. Yan and M. K. Wu Natl. Acad. Sci. U.S.A. 105 (2008) 14262.
  * [5] S. Margadonna, Y. Takabayashi, Y. Ohishi, Y. Mizuguchi, Y. Takano, T. Kagayama, T. Nakagawa, M. Takata, and K. Prassides Phys. Rev. B 80 (2009) 064506.
  * [6] M. Rotter, Phys. Rev. B  78 (2008) 020503(R).
  * [7] Z. Xu, J. Wen, G. Xu, Q. Jie, Z. Lin, Q. Li, S. Chi, D. K. Singh, G. Gu, and J. M. Tranquada Phys. Rev. B  82 (2010) 104525.
  * [8] C.Y. Moon, H. J. Choi, Phys. Rev. Lett.  104 (2010) 057003.
  * [9] F. Ma, W. Ji, J. Hu, Z. Y. Lu, T. Xiang Phys. Rev. Lett.  102 (2009) 177003.
  * [10] E. Z. Kuchinskii, I. A. Nekrasov, M. V. SadovskiiarXiv:1004.0801v1 (2010)
  * [11] and Yoshihiko Takano. J. Phys. Soc. Jpn.78 (2009) 074712.
  * [12] S. Li, C. De la Cruz, Q. Huang, Y. Chen, J. W. Lynn, J. Hu, Yi-Lin Huang, Fong-Chi Hsu, Kuo-Wei Yeh, Maw-Kuen Wu and P. Dai Phys. Rev. B  79 (2009) 054503.
  * [13] Y. Mizugucgi, Y. Hara, K. Deguchi, S. Tsuda, T. Yamaguchi, K. Takeda, H. Kotegawa, H. Tou and Y Takano Supercond. Sci. Technol.23 (2010) 054013.
  * [14] M. Bendele, P. Babkevich, S. Katrych, S. N. Gavasaliya, E. Pomjakushina, K. Conder, B. Roessli, A. T. Boothroyd, R. Khasanov and H. Keller Phys. Rev. B  82 (2010) 212504.
  * [15] R. Viennois, E. Giannini, D. van der Marel, R. cerny J. Solid State Chemistry 183 (2010) 769
  * [16] K. Nakayama, Y. Miyata, G. N. Phan, T Sato, Y. Tanabe, T. Urata, K. Tanigaki and T. Takahashi arXiv:1004.0857v1 (2014)
  * [17] Markus Aichhorn, Silke Biermann, Takashi Miyake, Antoine Georges and Masatoshi Imada Phys. Rev. B 82 (2010) 064504.
  * [18] Ansgar Liebsch Phys. Rev. B  84 (2011) 180505(R).
  * [19] P. Mishra, H. Lohani, R. A. Zargar, V. P. S. Awana and B. R. Sekhar J. Phys. Condens Matter 26 (2014) 425501.
  * [20] R. Yoshida, T. Wakita, H. Okazaki, Y. Mizuguchi, S. Tsuda, Y. Takano, H. Takeya, K. Hirata, T. Muro, M. Okawa, K. Ishizaka, S. Shin, H. Harima, M. Hirai, Y. Muraoka, and T. Yokoya J. Phys. Soc. Jpn. 78 (2009) 034708
  * [21] Kazumasa Horigane, Haruhiro Hiraka and Kenji Ohoyama J. Phys. Soc. Jpn. 78 (2009) 074718.
  * [22] Zargar R A, Pal A, Haﬁz A K, Awana V P S J. Supercond. Nov. Magn. 27 (2014) 897
  * [23] Takashi Miyake, Kazuma Nakamura, Ryotaro Arita and Masatoshi Imada J. Phys. Soc. Jpn. 79 (2010) 044705.
  * [24] O.K.Andersen, Phys. Rev. B 12 (1975) 3060.
  * [25] O.K.Andersen and O.Jepsen Phys. Rev. Lett. 53 (1984) 2571.
  * [26] R. W. Tank, O. Jepsen, A. Burkhardt and O. K. Andersen, The Stuttgart TB-LMTO-ASA programe,Version 47
  * [27] U. von Barth and L. Hedin, J. Phys. C, 4 (1971) 2064.
  * [28] Vladimir Anisimov, F. Aryasetiwan and A. I. Lichtenstein J. Phys. Condens Matter  9 (1997) 767.
  * [29] Giannozzi P et al. http://www.quantum-espresso.org.
  * [30] G. R. Grechnev, A. S. Panfilov, A. V. Fedorchenko, V.A. Desnenko, I.P. Zhuravleva, S.L. Gnatchenko, D.A. Chareev, O.S. Volkov, A.N. Vasiliev Ukr. J. Phys. 57 (2012) 171.
  * [31] L. Craco and S. Leoni EPL  92 (2010) 67003.
  * [32] Y. Xia, D. Qian, L. Wray, D. Hsieh, G.F. Chen, J.L. Luo, N.L. Wang and M.Z. Hasan Phys. Rev. Lett. 103 (2009) 037002.
  * [33] A. Yamaasaki, Y. Matsui, S.Imada, K. Takase, H. Azuma, T. Muro, Y.Kato, A.Higashiya, A.Sekiyama, S. Suga, M. Yabashi, K.Tamasaku, T. Ishikawa, K. Terashima, H. Kobori, A. Sugimura, N. Umeyama, H. Sato, Y. Hara, N. Miyagawa and S.I. Ikeda Phys. Rev. B  82 (2010) 184511.
  * [34] L. Craco and M.S. Laad arXiv:1001.3273v1, 2010
  * [35] A. Tamai, A.Y. Ganin, E. Rozbicki, J. Bacsa, W. Meevasana, P.D.C. King, M. Caffio, R. Schaub, S. Margadonna, K. Prassides, M.J. Rosseinsky and F. Baumberger Phys. Rev. Lett.  104 (2010) 097002.
  * [36] Alaska Subedi, Lijun Zhang, D.J. Singh and M.H. Du Phys. Rev. B 78 (2008) 134514.
  * [37] Fei Chen, Bo Zhou, Yan Zhang, Jia Wei Ou, Jia-Feng Zhao, Cheng He, Qing-Qin Ge, Masashi Arita, Kenya Shimada, Hirofumi Namatame, Masaki Taniguchi, Zhong-Yi Lu, Jiangping Hu, Xiao-Yu Cui and D.L. Feng Phys. Rev. B  81 (2010) 014526.
  * [38] L. Craco, M.S. Laad and S. Leoni arXiv:0910.3828v1, (2009)
  * [39] E. Ieki, K. Nakayama, Y. Miyata, T. Sato, H. Miao, N. Xu, X. -P. Wang, P. Zhang, T. Qian, P. Richard, Z. -J. Xu, J. S. Wen, G. D. Gu, H. Q. Luo, H. -H. Wen, H. Ding and T. Takahashi Phys. Rev. B  89 (2014) 140506(R).
  * [40] H. Miao, P. Richard, Y. Tanaka, K. Nakayama, T. Qian, K. Umezawa, T. Sato, Y.-M. Xu, Y. B. Shi, N. Xu, X.-P. Wang, P. Zhang, H.-B. Yang, Z.-J. Xu, J. S. Wen, G.-D. Gu, X. Dai, J.-P. Hu, T. Takahashi and H. Ding Phys. Rev. B  85 (2012) 094506.
  * [41] K. Nakayama, T. Sato, P. Richard, T. Kawahara, Y. Sekiba, T. Qian, G. F. Chen, J. L. Luo, N. L. Wang, H.Ding, and T. Takahashi Phys. Rev. Lett.  105 (2010) 197001.
  * [42] P. Starowicz, H. Schwab, J. Goraus, P. Zajdel, F. Forster, J. R. Rak, M. A. Green, I. Vobornik and F. Reinert arXiv:1304.7254v1, (2013)
  * [43] Antoine Georges, Luca De’ Medici and Jernej Mravlje arXiv:1207.3033v2, (2012)
  * [44] K Haule and G Kotliar New Journal of Physics  11 (2009) 025021
  * [45] Z. P. Yin, K. Haule and G. Kotliar Nature Materials Letters 10 (2011) 932
  * [46] Luca de’ Medici Phys. Rev. B  83 (2011) 205112.
  * [47] Nicola Lanata, Hugo U. R. Strand, Gianluca Giovannetti, Bo Hellsing, Luca de’ Medici and Massimo Capone Phys. Rev. B  87 (2013) 045122.
  * [48] Z. K. Liu, R.-H He, D.H. Lu, M. Yi, Y. L. Chen, M. Hashimoto, R. G. Moore, S.-K. Mo, E. A. Nowadnick, J. Hu, T. J. Liu, Z. Q. Mao, T. P. Devereaux, Z. Hussain and Z.-X. Shen Phys. Rev. Lett.  110 (2013) 037003.

Table 1:  In this table we have presented the change in the occupancy of
electrons in different d orbitals of Iron under the application of different
values of U and J in FeSe and FeTe compound. FeSe
---
U(eV) | J(eV) | dxy | dxz/yz | d${}_{3z^{2}-r^{2}}$ | d${}_{x^{2}-y^{2}}$
0.0 | 0.0 | 1.239 | 1.300 | 1.461 | 1.363
4.0 | 0.0 | 0.954 | 1.581 | 1.655 | 0.792
4.0 | 0.1 | 1.073 | 1.257 | 1.765 | 1.436
4.0 | 1.2 | 1.189 | 1.203 | 1.593 | 1.526
FeTe
0.0 | 0.0 | 1.242 | 1.308 | 1.462 | 1.363
3.5 | 0.0 | 1.182 | 1.109 | 1.569 | 1.637
3.5 | 0.1 | 1.174 | 1.138 | 1.558 | 1.634
3.5 | 1.2 | 1.189 | 1.209 | 1.514 | 1.572

Figure 1:  Total and partial DOS have been plotted for FeSe in (a) and for
FeTe in (1b). In (c) and (d) total and partial DOS have been calculated, where
Se-p and Te-p valence orbitals are doenfolded in FeSe and FeTe respectively.
Inset of (c) and (d) shows the anion tetrahedra in FeSe and FeTe respectively.
Fe, Se and Te atoms are denoted by red, blue and green clours respectively in
the inset pictures (c) and (d).

Figure 2:  Panel (a) shows changes in Fe-3d states in FeSe at different values
of U, where the value of U is written in red colour in each plot. Panel (b)
shows changes in Fe-3d states in FeTe at different values of U, where the
value of U is written in red colour in each plot.

Figure 3:  Panel (a) shows changes in Fe-3d states in FeSe at different values
of J, where the value of J is written in red colour in each plot. Panel (b)
shows changes in Fe-3d states in FeTe at different values of J, where the
value of J is written in red colour in each plot. The value of U is 4.0 and
3.5 eV in all FeSe and FeTe plots respectively.

Figure 4: Panel (a) shows changes in the DOS of different d orbitals of Fe
atom in FeSe at different values of U and J. Panel (B) shows changes in the
DOS of different d orbitals of Fe atom in FeTe at different values of U and J.
Different values of U and J are denoted by Green, Black, Red and Blue colours
at the top of the Figure.

Figure 5:  (a) and (b) shows band structure plot of FeSe at U = 0.0 and U =
4.0 eV respectively. (c) and (d) shows band structure plot of FeTe at U = 0.0
and U = 3.5 eV respectively. (e), (f), (g) and (h) shows fatness of bands
which are originated from dxy, dyz/xz, d${}_{x^{2}-y^{2}}$ and
d${}_{3z^{2}-r^{2}}$ respectively for FeSe at U and J = 0.0 eV.

Figure 6:  Panel (a) and (b) shows band structure plot for FeSe and FeTe at
different values of J = 0.0, 0.1, 0.8 and 1.2 eV from left to right
respectively. The value of U is 4.0 and 3.5 eV in all FeSe and FeTe plots
respectively.


# Fast Hankel Tensor-Vector Products and Application to Exponential Data
Fitting

Weiyang Ding E-mail<EMAIL_ADDRESS>School of Mathematical Sciences,
Fudan University, Shanghai, 200433, P. R. of China. This author is supported
by 973 Program Project under grant 2010CB327900. Liqun Qi Email:
<EMAIL_ADDRESS>Department of Applied Mathematics, The Hong Kong
Polytechnic University, Hung Hom, Kowloon, Hong Kong. This author is supported
by the Hong Kong Research Grant Council (Grant No. PolyU 502510, 502111,
501212, 501913). Yimin Wei E-mail<EMAIL_ADDRESS>and<EMAIL_ADDRESS>School of Mathematical Sciences and Key Laboratory of Mathematics for
Nonlinear Sciences, Fudan University, Shanghai, 200433, P. R. of China. This
author is supported by the National Natural Science Foundation of China under
grant 11271084.

###### Abstract

This paper is contributed to a fast algorithm for Hankel tensor-vector
products. For this purpose, we first discuss a special class of Hankel tensors
that can be diagonalized by the Fourier matrix, which is called _anti-
circulant_ tensors. Then we obtain a fast algorithm for Hankel tensor-vector
products by embedding a Hankel tensor into a larger anti-circulant tensor. The
computational complexity is about $\mathcal{O}(m^{2}n\log mn)$ for a square
Hankel tensor of order $m$ and dimension $n$, and the numerical examples also
show the efficiency of this scheme. Moreover, the block version for multi-
level block Hankel tensors is discussed as well. Finally, we apply the fast
algorithm to exponential data fitting and the block version to 2D exponential
data fitting for higher performance.

Keywords: Hankel tensor, block Hankel tensor, anti-circulant tensor, fast
tensor-vector product, fast Fourier transformation, higher-order singular
value decomposition, exponential data fitting.

2000 AMS Subject Classification: 15A69, 15B05, 65T50, 68M07.

##  1 Introduction

Hankel structures arise frequently in signal processing [17]. Besides Hankel
matrices, tensors with different Hankel structures have also been applied. As
far as we know, the concept of “Hankel tensor” was first introduced in [15] by
Luque and Thibon. Papy and Boyer et al. employed Hankel tensors and tensors
with other Hankel structures in exponential data fitting in [2, 18, 19]. Then
Badeau and Boyer discussed the higher-order singular value decompositions
(HOSVD) of some structured tensors, including Hankel tensors, in [1].
Recently, Qi investigated the spectral properties of Hankel tensor in [23]
largely via the generating function, Song and Qi proposed some spectral
properties of Hilbert tensors in [25], which are special Hankel tensors.

A tensor $\mathcal{H}\in\mathbb{C}^{n_{1}\times n_{2}\times\dots\times n_{m}}$
is called a _Hankel tensor_ , if

$\mathcal{H}_{i_{1}i_{2}\dots i_{m}}=\phi(i_{1}+i_{2}+\dots+i_{m})$

for all $i_{k}=0,1,\dots,n_{k}-1$ ($k=1,2,\dots,m$). We call $\mathcal{H}$ a
_square Hankel tensor_ , when $n_{1}=n_{2}=\dots=n_{m}$. Note that the degree
of freedom of a Hankel tensor is
$d_{\mathcal{H}}:=n_{1}+n_{2}+\dots+n_{m}-m+1$. Thus a vector ${\bf h}$ of
length $d_{\mathcal{H}}$, which is called the _generating vector_ of
$\mathcal{H}$, defined by

$h_{k}=\phi(k),\ k=0,1,\dots,d_{\mathcal{H}}-1$

can totally determine the Hankel tensor $\mathcal{H}$, when the tensor size is
fixed. Further, if the entries of ${\bf h}$ can be written as

$h_{k}=\int_{-\infty}^{+\infty}t^{k}f(t){\rm d}t,$

we then call $f(t)$ the _generating function_ of $\mathcal{H}$. The generating
function of a square Hankel tensor is essential for studying eigenvalues,
positive semi-definiteness, and copositiveness etc. (cf. [23]).

The fast algorithm for Hankel or Toeplitz matrix-vector products using fast
Fourier transformations (FFT) is well-known (cf. [4, 16]). However the topics
on Hankel tensor computations are seldom discussed. We propose an analogous
fast algorithm for Hankel tensor-vector products and application to signal
processing. In Section 2, we begin with a special class of Hankel tensors
called _anti-circulant_ tensors. By using some properties of anti-circulant
tensors, the Hankel tensor-vector products can be performed in complexity
$\mathcal{O}(md_{\mathcal{H}}\log d_{\mathcal{H}})$, which is described in
Section 3. Then in Section 4, we apply the fast algorithm and its block
version to 1D exponential data fitting and 2D exponential data fitting.
Finally, we draw some conclusions and display some future research in Section
5.

##  2 Anti-Circulant Tensor

Circulant matrix is famous, which is a special class of Toeplitz matrices [4,
16]. The first column entries of a circulant matrix shift down when moving
right, as shown in the following $3$-by-$3$ example

$\begin{bmatrix}c_{0}&c_{2}&c_{1}\\\ c_{1}&c_{0}&c_{2}\\\
c_{2}&c_{1}&c_{0}\end{bmatrix}.$

If the first column entries of a matrix shift up when moving right, as shown
in the following

$\begin{bmatrix}c_{0}&c_{1}&c_{2}\\\ c_{1}&c_{2}&c_{0}\\\
c_{2}&c_{0}&c_{1}\end{bmatrix},$

then it is a special Hankel matrix. And we name it an _anti-circulant matrix_.
Naturally, we generalize the anti-circulant matrix to the tensor case. A
square Hankel tensor $\mathcal{C}$ of order $m$ and dimension $n$ is called an
_anti-circulant tensor_ , if its generating vector ${\bf h}$ satisfies that

$h_{k}=h_{l},\quad\text{if }k\equiv l\,({\rm mod}\,n).$

Thus the generating vector is periodic and displayed as

${\bf h}=(\underbrace{h_{0},h_{1},\dots,h_{n-1}}_{{\bf
c}^{\top}},\underbrace{h_{n},h_{n+1},\dots,h_{2n-1}}_{{\bf
c}^{\top}},\dots,\underbrace{h_{(m-1)n},\dots,h_{m(n-1)}}_{{\bf
c}(0:n-m)^{\top}})^{\top}.$

Since the vector ${\bf c}$, which is exactly the “first” column
$\mathcal{C}(:,0,\cdots,0)$, contains all the information about $\mathcal{C}$
and is more compact than the generating vector, we call it the _compressed
generating vector_ of the anti-circulant tensor. For instance, a $3\times
3\times 3$ anti-circulant tensor $\mathcal{C}$ is unfolded by mode-$1$ into

${\rm Unfold}_{1}(\mathcal{C})=\left[\begin{array}[]{ccc}c_{0}&c_{1}&c_{2}\\\
c_{1}&c_{2}&c_{0}\\\
c_{2}&c_{0}&c_{1}\end{array}\right.\left|\begin{array}[]{ccc}c_{1}&c_{2}&c_{0}\\\
c_{2}&c_{0}&c_{1}\\\
c_{0}&c_{1}&c_{2}\end{array}\right|\left.\begin{array}[]{ccc}c_{2}&c_{0}&c_{1}\\\
c_{0}&c_{1}&c_{2}\\\ c_{1}&c_{2}&c_{0}\end{array}\right],$

and its compressed generating vector is ${\bf c}=[c_{0},c_{1},c_{2}]^{\top}$.
Note that the degree of freedom of an anti-circulant tensor is always $n$ no
matter how large its order $m$ will be.

###  2.1 Diagonalization

One of the essential properties of circulant matrices is that every circulant
matrix can be diagonalized by the Fourier matrix, where the Fourier matrix of
size $n$ is defined as $F_{n}=\big{(}\exp(-\frac{2\pi{\bf
i}}{n}jk)\big{)}_{j,k=0,1,\dots,n-1}$. Actually, the Fourier matrix is exactly
the Vandermonde matrix for the roots of unity, and it is also a unitary matrix
up to the normalization factor

$F_{n}F_{n}^{\ast}=F_{n}^{\ast}F_{n}=nI_{n},$

where $I_{n}$ is the identity matrix of $n\times n$ and $F_{n}^{\ast}$ is the
conjugate transpose of $F_{n}$. We will show that anti-circulant tensors also
have a similar property, which brings much convenience for both analysis and
computations.

In order to describe this property, we recall the definition of mode-$p$
tensor-matrix product first. It should be pointed out that the tensor-matrix
products in this paper are slightly different with some standard notations
(cf. [6, 12]) just for easy use and simple descriptions. Let $\mathcal{A}$ be
a tensor of $n_{1}\times n_{2}\times\dots\times n_{m}$ and $M$ be a matrix of
$n_{p}\times n_{p}^{\prime}$, then the mode-$p$ product
$\mathcal{A}\times_{p}M$ is another tensor of $n_{1}\times\dots\times
n_{p-1}\times n_{p}^{\prime}\times n_{p+1}\times\dots\times n_{m}$ defined as

$(\mathcal{A}\times_{p}M)_{i_{1}\dots i_{p-1}ji_{p+1}\dots
i_{m}}=\sum_{i_{p}=0}^{n_{p}-1}\mathcal{A}_{i_{1}\dots
i_{p-1}i_{p}i_{p+1}\dots i_{m}}M_{i_{p}j}.$

In the standard notation system, two indices in “$M_{i_{p}j}$” should be
exchanged. There are some basic properties of the tensor-matrix products

  1. 1.

$\mathcal{A}\times_{p}M_{p}\times_{q}M_{q}=\mathcal{A}\times_{q}M_{q}\times_{p}M_{p}$,
if $p\neq q$,

  2. 2.

$\mathcal{A}\times_{p}M_{p_{1}}\times_{p}M_{p_{2}}=\mathcal{A}\times_{p}(M_{p_{1}}M_{p_{2}})$,

  3. 3.

$\mathcal{A}\times_{p}M_{p_{1}}+\mathcal{A}\times_{p}M_{p_{2}}=\mathcal{A}\times_{p}(M_{p_{1}}+M_{p_{2}})$,

  4. 4.

$\mathcal{A}_{1}\times_{p}M+\mathcal{A}_{2}\times_{p}M=(\mathcal{A}_{1}+\mathcal{A}_{2})\times_{p}M$.

Particularly, when $\mathcal{A}$ is a matrix, the mode-$1$ and mode-$2$
products can be written as

$\mathcal{A}\times_{1}M_{1}\times_{2}M_{2}=M_{1}^{\top}\mathcal{A}M_{2}.$

Notice that this is totally different with $M_{1}^{\ast}\mathcal{A}M_{2}$!
($M_{1}^{\ast}$ is the conjugate transpose of $M_{1}$.) We will adopt some
notations from Qi’s paper [23]

$\begin{split}\mathcal{A}{\bf x}^{m-1}&=\mathcal{A}\hskip 26.0pt\times_{2}{\bf
x}\dots\times_{m}{\bf x},\\\ \mathcal{A}{\bf x}^{m}&=\mathcal{A}\times_{1}{\bf
x}\times_{2}{\bf x}\dots\times_{m}{\bf x}.\end{split}$

We are now ready to state our main result about anti-circulant tensors.

######  Theorem 2.1.

A square tensor of order $m$ and dimension $n$ is an anti-circulant tensor if
and only if it can be diagonalized by the Fourier matrix of size $n$, that is,

$\mathcal{C}=\mathcal{D}F_{n}^{m}:=\mathcal{D}\times_{1}F_{n}\times_{2}F_{n}\dots\times_{m}F_{n},$

where $\mathcal{D}$ is a diagonal tensor and ${\rm diag}(\mathcal{D})={\rm
ifft}({\bf c})$. Here, “${\rm ifft}$” is a Matlab-type symbol, an abbreviation
of inverse fast Fourier transformation.

###### Proof.

It is direct to verify that a tensor that can be expressed as
$\mathcal{D}F_{n}^{m}$ is anti-circulant. Thus we only need to prove that
every anti-circulant tensor can be written into this form. And this can be
done constructively.

First, assume that an anti-circulant tensor $\mathcal{C}$ could be written
into $\mathcal{D}F_{n}^{m}$. Then how do we obtain the diagonal entries of
$\mathcal{D}$ from $\mathcal{C}$? Since

$\begin{split}{\rm diag}(\mathcal{D})&=\mathcal{D}{\bf
1}^{m-1}=\frac{1}{n^{m}}\big{(}\mathcal{C}(F_{n}^{\ast})^{m}\big{)}{\bf
1}^{m-1}=\frac{1}{n^{m}}\bar{F}_{n}\big{(}\mathcal{C}(F_{n}^{\ast}{\bf
1})^{m-1}\big{)}\\\ &=\frac{1}{n}\bar{F}_{n}(\mathcal{C}{\bf
e}_{0}^{m-1})=\frac{1}{n}\bar{F}_{n}{\bf c},\end{split}$

where ${\bf 1}=[1,1,\dots,1]^{\top}$, ${\bf e}_{0}=[1,0,\dots,0]^{\top}$, and
${\bf c}$ is the compressed generating vector of $\mathcal{C}$, then the
diagonal entries of $\mathcal{D}$ can be computed by an inverse fast Fourier
transformation (IFFT)

${\rm diag}(\mathcal{D})={\rm ifft}({\bf c}).$

Finally, it is enough to check that $\mathcal{C}=\mathcal{D}F^{m}$ with ${\rm
diag}(\mathcal{D})={\rm ifft}({\bf c})$ directly. Therefore, every anti-
circulant tensor is diagonalized by the Fourier matrix of proper size. ∎

From the expression $\mathcal{C}=\mathcal{D}F_{n}^{m}$, we have a corollary
about the spectra of anti-circulant tensors. The definitions of tensor
Z-eigenvalues and H-eigenvalues follow the ones in [21, 22].

######  Corollary 2.2.

An anti-circulant tensor $\mathcal{C}$ of order $m$ and dimension $n$ with the
compressed generating vector ${\bf c}$ has a Z-eigenvalue/H-eigenvalue
$n^{m-2}{\bf 1}^{\top}{\bf c}$ with eigenvector ${\bf 1}$. When $n$ is even,
it has another Z-eigenvalue $n^{m-2}\tilde{\bf 1}^{\top}{\bf c}$ with
eigenvector $\tilde{\bf 1}$, where $\tilde{\bf 1}=[1,-1,\dots,1,-1]^{\top}$;
Moreover, this is also an H-eigenvalue if $m$ is even.

###### Proof.

Check it directly:

$\begin{split}\mathcal{C}{\bf
1}^{m-1}&=F_{n}^{\top}\big{(}\mathcal{D}(F_{n}{\bf
1})^{m-1}\big{)}=n^{m-1}F_{n}^{\top}(\mathcal{D}{\bf
e}_{0}^{m-1})=n^{m-1}\mathcal{D}_{1,1,\dots,1}\cdot F_{n}^{\top}{\bf e}_{0}\\\
&=n^{m-2}({\bf e}_{0}^{\top}\bar{F}_{n}{\bf c}){\bf 1}=(n^{m-2}{\bf
1}^{\top}{\bf c}){\bf 1}.\end{split}$

The proof of the rest part is similar, so we omit it. ∎

###  2.2 Block Tensors

Block structures arise in a variety of applications in scientific computing
and engineering (cf. [11, 17]). And block matrices are successfully applied to
many disciplines. Surely, we can also define block tensors. A _block tensor_
is understood as a tensor whose entries are also tensors. Let $\mathcal{A}$ be
a block tensor. Then the $(j_{1},j_{2},\dots,j_{m})$-th block is denoted as
$\mathcal{A}^{(j_{1}j_{2}\dots j_{m})}$, and the
$(i_{1},i_{2},\dots,i_{m})$-th entry of the $(j_{1},j_{2},\dots,j_{m})$-th
block is denoted as $\mathcal{A}_{i_{1}i_{2}\dots i_{m}}^{(j_{1}j_{2}\dots
j_{m})}$.

If a block tensor can be regarded as a Hankel tensor or an anti-circulant
tensor with tensor entries, then we call it a _block Hankel tensor_ or a
_block anti-circulant tensor_ , respectively. Moreover, its generating vector
${\bf h}^{(b)}$ or compressed generating vector ${\bf c}^{(b)}$ with tensor
entries is called the _block generating vector_ or _block compressed
generating vector_ , respectively. For instance, let $\mathcal{H}$ be a block
Hankel tensor of size
$(n_{1}N_{1})\times(n_{2}N_{2})\times\dots\times(n_{m}N_{m})$ with blocks of
size $n_{1}\times n_{2}\times\dots\times n_{m}$, then ${\bf h}^{(b)}$ is a
block tensor of size $[(N_{1}+N_{2}+\cdots+N_{m}-m+1)n_{1}]\times
n_{2}\times\dots\times n_{m}$. Recall the definition of Kronecker product [10]

$A\otimes B=\begin{bmatrix}A_{11}B&A_{12}B&\cdots&A_{1q}B\\\
\vdots&\vdots&\ddots&\vdots\\\ A_{p1}B&A_{p2}B&\cdots&A_{pq}B\end{bmatrix},$

where $A$ and $B$ are two matrices of arbitrary sizes. Then it can be proved
following Theorem 2.1 that a block anti-circulant tensor $\mathcal{C}$ can be
block-diagonalized by $F_{N}\otimes I$, that is,

$\mathcal{C}=\mathcal{D}^{(b)}(F_{N}\otimes I)^{m},$

where $\mathcal{D}^{(b)}$ is a block diagonal tensor with diagonal blocks
${\bf c}^{(b)}\times_{1}\big{(}\frac{1}{N}\bar{F}_{N}\otimes I\big{)}$ and
$\bar{F}_{N}$ is the conjugate of $F_{N}$.

Furthermore, when the blocks of a block Hankel tensor are also Hankel tensors,
we call it a _block Hankel tensor with Hankel blocks_ , or BHHB tensor for
short. Then its block generating vector can be reduced to a matrix, which is
called the _generating matrix_ $H$ of a BHHB tensor

$H=[{\bf h}_{0},{\bf h}_{1},\dots,{\bf
h}_{N_{1}+\dots+N_{m}-m}]\in\mathbb{C}^{(n_{1}+n_{2}+\dots+n_{m}-m+1)\times(N_{1}+N_{2}+\dots+N_{m}-m+1)},$

where ${\bf h}_{k}$ is the generating vector of the $k$-th Hankel block in
${\bf h}^{(b)}$. Similarly, when the blocks of a block anti-circulant tensor
are also anti-circulant tensors, we call it a _block anti-circulant tensor
with anti-circulant blocks_ , or BAAB tensor for short. Then its _compressed
generating matrix_ $C$ is defined as

$C=[{\bf c}_{0},{\bf c}_{1},\dots,{\bf c}_{N-1}]\in\mathbb{C}^{n\times N},$

where ${\bf c}_{k}$ is the compressed generating vector of the $k$-th anti-
circulant block in the block compressed generating vector ${\bf c}^{(b)}$. We
can also verify that a BAAB tensor $\mathcal{C}$ can be diagonalized by
$F_{N}\otimes F_{n}$, that is,

$\mathcal{C}=\mathcal{D}(F_{N}\otimes F_{n})^{m},$

where $\mathcal{D}$ is a diagonal tensor with diagonal ${\rm
diag}(\mathcal{D})=\frac{1}{nN}{\rm vec}(\bar{F}_{n}C\bar{F}_{N})$, which can
be computed by 2D inverse fast Fourier transformation (IFFT2). Here ${\rm
vec}(\cdot)$ denotes the vectorization operator (cf. [10]).

We can even define higher-level block Hankel tensors. For instance, a block
Hankel tensor with BHHB blocks is called a _level- $3$ block Hankel tensor_,
and it is easily understood that a level-$3$ block Hankel tensor has the
generating tensor of order $3$. Generally, a block Hankel or anti-circulant
tensor with level-$(k$-$1)$ block Hankel or anti-circulant blocks is called a
_level- $k$ block Hankel or anti-circulant tensor_, respectively. Furthermore,
a level-$k$ block anti-circulant tensor $\mathcal{C}$ can be diagonalized by
$F_{n^{(k)}}\otimes F_{n^{(k-1)}}\otimes\dots\otimes F_{n^{(1)}}$, that is,

$\mathcal{C}=\mathcal{D}(F_{n^{(k)}}\otimes F_{n^{(k-1)}}\otimes\dots\otimes
F_{n^{(1)}})^{m},$

where $\mathcal{D}$ is a diagonal tensor with diagonal that can be computed by
multi-dimensional inverse fast Fourier transformation.

##  3 Fast Hankel Tensor-Vector Product

General tensor-vector products without structures are very expensive for the
high order and the large size that a tensor could be of. For a square tensor
$\mathcal{A}$ of order $m$ and dimension $n$, the computational complexity of
a tensor-vector product $\mathcal{A}{\bf x}^{m-1}$ or $\mathcal{A}{\bf x}^{m}$
is $\mathcal{O}(n^{m})$. However, since Hankel tensors and anti-circulant
tensors have very low degrees of freedom, it can be expected that there is a
much faster algorithm for Hankel tensor-vector products. We focus on the
following two types of tensor-vector products

${\bf y}=\mathcal{A}\times_{2}{\bf x}_{2}\cdots\times_{m}{\bf x}_{m}\text{ and
}\alpha=\mathcal{A}\times_{1}{\bf x}_{1}\times_{2}{\bf
x}_{2}\cdots\times_{m}{\bf x}_{m},$

which will be extremely useful to applications.

The fast algorithm for anti-circulant tensor-vector products is easy to derive
from Theorem 2.1. Let $\mathcal{C}=\mathcal{D}F_{n}^{m}$ be an anti-circulant
tensor of order $m$ and dimension $n$ with the compressed generating vector
${\bf c}$. Then for vectors ${\bf x}_{2},{\bf x}_{3},\dots,{\bf
x}_{m}\in\mathbb{C}^{n}$, we have

${\bf y}=\mathcal{C}\times_{2}{\bf x}_{2}\cdots\times_{m}{\bf
x}_{m}=F_{n}(\mathcal{D}\times_{2}F_{n}{\bf x}_{2}\cdots\times_{m}F_{n}{\bf
x}_{m}).$

Recall that ${\rm diag}(\mathcal{D})={\rm ifft}({\bf c})$ and $F_{n}{\bf
v}={\rm fft}({\bf v})$, where “${\rm fft}$” is a Matlab-type symbol, an
abbreviation of fast Fourier transformation. So the fast procedure for
computing the vector ${\bf y}$ is

${\bf y}={\rm fft}\big{(}{\rm ifft}({\bf c}).\ast{\rm fft}({\bf
x}_{2}).\ast\cdots.\ast{\rm fft}({\bf x}_{m})\big{)},$

where ${\bf u}.\ast{\bf v}$ multiplies two vectors element-by-element.
Similarly, for vectors ${\bf x}_{1},{\bf x}_{2},\dots,{\bf
x}_{m}\in\mathbb{C}^{n}$, we have

$\alpha=\mathcal{C}\times_{1}{\bf x}_{1}\times_{2}{\bf
x}_{2}\cdots\times_{m}{\bf x}_{m}=\mathcal{D}\times_{1}F_{n}{\bf
x}_{1}\times_{2}F_{n}{\bf x}_{2}\cdots\times_{m}F_{n}{\bf x}_{m},$

and the fast procedure for computing the scalar $\alpha$ is

$\alpha={\rm ifft}({\bf c})^{\top}\big{(}{\rm fft}({\bf x}_{1}).\ast{\rm
fft}({\bf x}_{2}).\ast\cdots.\ast{\rm fft}({\bf x}_{m})\big{)}.$

Since the computational complexity of either FFT or IFFT is $\mathcal{O}(n\log
n)$, both two types of anti-circulant tensor-vector products can be obtained
with complexity $\mathcal{O}\big{(}(m+1)n\log n\big{)}$, which is much faster
than the product of a general $n$-by-$n$ matrix with a vector.

For deriving the fast algorithm for Hankel tensor-vector products, we embed a
Hankel tensor into a larger anti-circulant tensor. Let
$\mathcal{H}\in\mathbb{C}^{n_{1}\times n_{2}\times\dots\times n_{m}}$ be a
Hankel tensor with the generating vector ${\bf h}$. Denote
$\mathcal{C}_{\mathcal{H}}$ as the anti-circulant tensor of order $m$ and
dimension $d_{\mathcal{H}}=n_{1}+n_{2}+\dots+n_{m}-m+1$ with the compressed
generating vector ${\bf h}$. Then we will find out that $\mathcal{H}$ is in
the “upper left frontal” corner of $\mathcal{C}_{\mathcal{H}}$ as shown in
Figure 1.

Figure 1: Embed a Hankel tensor into an anti-circulant tensor.

Hence, we have

$\begin{split}&\mathcal{C}_{\mathcal{H}}\times_{2}\begin{bmatrix}{\bf
x}_{2}\\\ {\bf 0}\end{bmatrix}\dots\times_{m}\begin{bmatrix}{\bf x}_{m}\\\
{\bf 0}\end{bmatrix}=\begin{bmatrix}\mathcal{H}\times_{2}{\bf
x}_{2}\dots\times_{m}{\bf x}_{m}\\\ {\bf{\dagger}}\end{bmatrix},\\\
&\mathcal{C}_{\mathcal{H}}\times_{1}\begin{bmatrix}{\bf x}_{1}\\\ {\bf
0}\end{bmatrix}\dots\times_{m}\begin{bmatrix}{\bf x}_{m}\\\ {\bf
0}\end{bmatrix}=\mathcal{H}\times_{1}{\bf x}_{1}\dots\times_{m}{\bf
x}_{m},\end{split}$

so that the Hankel tensor-vector products can be realized by multiplying a
larger anti-circulant tensor by some augmented vectors. Therefore, the fast
procedure for computing ${\bf y}=\mathcal{H}\times_{2}{\bf
x}_{2}\dots\times_{m}{\bf x}_{m}$ is

$\left\\{\begin{array}[]{l}\tilde{\bf x}_{p}=[{\bf
x}_{p}^{\top},\underbrace{0,0,\dots,0}_{d_{\mathcal{H}}-n_{p}}]^{\top},\
p=2,3,\dots,m,\\\ \tilde{\bf y}={\rm fft}\big{(}{\rm ifft}({\bf h}).\ast{\rm
fft}(\tilde{\bf x}_{2}).\ast\cdots.\ast{\rm fft}(\tilde{\bf x}_{m})\big{)},\\\
{\bf y}=\tilde{\bf y}(0:n_{1}-1),\end{array}\right.$

and the fast procedure for computing $\alpha=\mathcal{H}\times_{1}{\bf
x}_{1}\times_{2}{\bf x}_{2}\dots\times_{m}{\bf x}_{m}$ is

$\left\\{\begin{array}[]{l}\tilde{\bf x}_{p}=[{\bf
x}_{p}^{\top},\underbrace{0,0,\dots,0}_{d_{\mathcal{H}}-n_{p}}]^{\top},\
p=1,2,\dots,m,\\\ \alpha={\rm ifft}({\bf h})^{\top}\big{(}{\rm fft}(\tilde{\bf
x}_{1}).\ast{\rm fft}(\tilde{\bf x}_{2}).\ast\cdots.\ast{\rm fft}(\tilde{\bf
x}_{m})\big{)}.\end{array}\right.$

Moreover, the computational complexity is
$\mathcal{O}\big{(}(m+1)d_{\mathcal{H}}\log d_{\mathcal{H}}\big{)}$. When the
Hankel tensor is a square tensor, the complexity is at the level
$\mathcal{O}(m^{2}n\log mn)$, which is much smaller than the complexity
$\mathcal{O}(n^{m})$ of non-structured products.

Apart from the low computational complexity, our algorithm for Hankel tensor-
vector products has two advantages. One is that this scheme is compact, that
is, there is no redundant element in the procedure. It is not required to form
the Hankel tensor explicitly. Just the generating vector is needed. Another
advantage is that our algorithm treats the tensor as an ensemble instead of
multiplying the tensor by vectors mode by mode.

######  Example 3.1.

We choose $3^{\rm th}$-order square Hankel tensors $\mathcal{H}$ of different
sizes. We compute the tensor-vector products $\mathcal{H}{\bf x}^{m-1}$ by
using our fast algorithm and the non-structured algorithm directly based on
the definition. The comparison of the running times is shown in Figure 2. From
the results, we can see that the running time of our algorithm increases far
slowly than that of the non-structured algorithm just as the theoretical
analysis. Moreover, the difference in running times is not only the low
computational complexity, but also the absence of forming the Hankel tensor
explicitly in our algorithm.

Figure 2: The running times of 100 products for each size and algorithm.

For BAAB and BHHB cases, we also have fast algorithms for the tensor-vector
products. Let $\mathcal{C}$ be a BAAB tensor of order $m$ with the compressed
generating matrix $C\in\mathbb{C}^{n\times N}$. Since $\mathcal{C}$ can be
diagonalized by $F_{N}\otimes F_{n}$, i.e.,

$\mathcal{C}=\mathcal{D}(F_{N}\otimes F_{n})^{m},$

we have for vectors ${\bf x}_{2},{\bf x}_{3},\dots,{\bf
x}_{m}\in\mathbb{C}^{nN}$

${\bf y}=\mathcal{C}\times_{2}{\bf x}_{2}\cdots\times_{m}{\bf
x}_{m}=(F_{N}\otimes F_{n})\big{(}\mathcal{D}\times_{2}(F_{N}\otimes
F_{n}){\bf x}_{2}\cdots\times_{m}(F_{N}\otimes F_{n}){\bf x}_{m}\big{)}.$

Recall the vectorization operator and its inverse operator

$\begin{split}{\rm
vec}(A)&=[A_{:,0}^{\top},A_{:,1}^{\top},\dots,A_{:,N-1}^{\top}]^{\top}\in\mathbb{C}^{nN},\\\
{\rm vec}_{n,N}^{-1}({\bf v})&=[{\bf v}_{0:n-1},{\bf v}_{n:2n-1},\dots,{\bf
v}_{(N-1)n:Nn-1}]\in\mathbb{C}^{n\times N},\end{split}$

for matrix $A\in\mathbb{C}^{n\times N}$ and vector ${\bf
v}\in\mathbb{C}^{nN}$, and the relation holds

$(B\otimes A){\bf v}={\rm vec}\big{(}A\cdot{\rm vec}_{n,N}^{-1}({\bf v})\cdot
B^{\top}\big{)}.$

So $(F_{N}\otimes F_{n}){\bf x}_{p}={\rm vec}\big{(}F_{n}\cdot{\rm
vec}_{n,N}^{-1}({\bf x}_{p})\cdot F_{N}\big{)}$ can be computed by 2D fast
Fourier transformation (FFT2). Then the fast procedure for computing ${\bf
y}=\mathcal{C}\times_{2}{\bf x}_{2}\dots\times_{m}{\bf x}_{m}$ is

$\left\\{\begin{array}[]{l}X_{p}={\rm vec}_{n,N}^{-1}({\bf x}_{p}),\
p=2,3,\dots,m,\\\ Y={\rm fft2}\big{(}{\rm ifft2}(C).\ast{\rm
fft2}(X_{2}).\ast\cdots.\ast{\rm fft2}(X_{m})\big{)},\\\ {\bf y}={\rm
vec}(Y),\end{array}\right.$

and the fast procedure for computing $\alpha=\mathcal{C}\times_{1}{\bf
x}_{1}\times_{2}{\bf x}_{2}\dots\times_{m}{\bf x}_{m}$ is

$\left\\{\begin{array}[]{l}X_{p}={\rm vec}_{n,N}^{-1}({\bf x}_{p}),\
p=1,2,\dots,m,\\\ \alpha=\big{\langle}{\rm ifft2}(C),{\rm
fft2}(X_{1}).\ast{\rm fft2}(X_{2}).\ast\cdots.\ast{\rm
fft2}(X_{m})\big{\rangle},\end{array}\right.$

where $\langle A,B\rangle$ denotes

$\langle A,B\rangle=\sum_{j,k}A_{jk}B_{jk}.$

For a BHHB tensor $\mathcal{H}$ with the generating matrix $H$, we do the
embedding twice. First we embed each Hankel block into a larger anti-circulant
block, and then we embed the block Hankel tensor with anti-circulant blocks
into a BAAB tensor $\mathcal{C}_{\mathcal{H}}$. Notice that the compressed
generating matrix of $\mathcal{C}_{\mathcal{H}}$ is exactly the generating
matrix of $\mathcal{H}$. Hence we have the fast procedure for computing ${\bf
y}=\mathcal{H}\times_{2}{\bf x}_{2}\dots\times_{m}{\bf x}_{m}$

$\left\\{\begin{array}[]{l}\tilde{X}_{p}=\underbrace{\begin{bmatrix}{\rm
vec}_{n_{p},N_{p}}^{-1}({\bf x}_{p})&{\rm O}\\\ {\rm O}&{\rm
O}\end{bmatrix}}_{N_{1}+N_{2}+\cdots+N_{m}-m+1}\left.\begin{matrix}\ \\\ \
\end{matrix}\right\\}{}^{n_{1}+n_{2}+\cdots+n_{m}-m+1},\ p=2,3,\dots,m,\\\
\tilde{Y}={\rm fft2}\big{(}{\rm ifft2}(H).\ast{\rm
fft2}(\tilde{X}_{2}).\ast\cdots.\ast{\rm fft2}(\tilde{X}_{m})\big{)},\\\ {\bf
y}={\rm vec}\big{(}\tilde{Y}(0:n_{1}-1,0:N_{1}-1)\big{)}.\end{array}\right.$

Sometimes in applications there is no need to do the vectorization in the last
line, and we just keep it as a matrix for later use. We also have the fast
procedure for computing $\alpha=\mathcal{H}\times_{1}{\bf x}_{1}\times_{2}{\bf
x}_{2}\dots\times_{m}{\bf x}_{m}$

$\left\\{\begin{array}[]{l}\tilde{X}_{p}=\underbrace{\begin{bmatrix}{\rm
vec}_{n_{p},N_{p}}^{-1}({\bf x}_{p})&{\rm O}\\\ {\rm O}&{\rm
O}\end{bmatrix}}_{N_{1}+N_{2}+\cdots+N_{m}-m+1}\left.\begin{matrix}\ \\\ \
\end{matrix}\right\\}{}^{n_{1}+n_{2}+\cdots+n_{m}-m+1},\ p=2,3,\dots,m,\\\
\alpha=\big{\langle}{\rm ifft2}(H),{\rm fft2}(\tilde{X}_{1}).\ast{\rm
fft2}(\tilde{X}_{2}).\ast\cdots.\ast{\rm
fft2}(\tilde{X}_{m})\big{\rangle}.\end{array}\right.$

Similarly, we can also derive the fast algorithms for higher-level block
Hankel tensor-vector products using the multi-dimensional FFT.

##  4 Exponential Data Fitting

Exponential data fitting is very important in many applications in scientific
computing and engineering, which represents the signals as a sum of
exponentially damped sinusoids. The computations and applications of
exponential data fitting are generally studied, and the reader is interested
in these topics can refer to [20]. Furthermore, Papy et al. introduced a
higher-order tensor approach into exponential data fitting in [18, 19] by
connecting it with the Vandermonde decomposition of a Hankel tensor. And they
showed that the tensor approach can do better than the classical matrix
approach.

###  4.1 Hankel Tensor Approach

Assume that we get a one-dimensional noiseless signal with $N$ complex samples
$\\{x_{n}\\}_{n=0}^{N-1}$, and this signal is modelled as a sum of $K$
exponentially damped complex sinusoids, i.e.,

$x_{n}=\sum_{k=1}^{K}a_{k}\exp({\bf i}\varphi_{k})\exp\big{(}(-\alpha_{k}+{\bf
i}\omega_{k})n\Delta t\big{)},$

where ${\bf i}=\sqrt{-1}$, $\Delta t$ is the sampling interval, and the
amplitudes $a_{k}$, the phases $\varphi_{k}$, the damping factors
$\alpha_{k}$, and the pulsations $\omega_{k}$ are the parameters of the model.
The signal can also be expressed as

$x_{n}=\sum_{k=1}^{K}c_{k}z_{k}^{n},$

where $c_{k}=a_{k}\exp({\bf i}\varphi_{k})$ and
$z_{k}=\exp\big{(}(-\alpha_{k}+{\bf i}\omega_{k})\Delta t\big{)}$. Here
$c_{k}$ is called the $k$-th complex amplitude including the phase, and
$z_{k}$ is called the $k$-th pole of the signal. A part of the aim of
exponential data fitting is to estimate the poles $\\{z_{k}\\}_{k=1}^{K}$ from
the data $\\{x_{n}\\}_{n=0}^{N-1}$.

Denote vector ${\bf x}=(x_{n})_{n=0}^{N-1}$. Then we construct a Hankel tensor
$\mathcal{H}$ of a fixed order $m$ and size $I_{1}\times
I_{2}\times\dots\times I_{m}$ with the generating vector ${\bf x}$. The order
$m$ can be chosen arbitrarily and the size $I_{p}$ of each dimension should be
no less than $K$. And we always deal with $3$-order tensors in this paper.
Moreover, it is advised to construct a Hankel tensor as close to a square
tensor as well (cf. [18, 19, 20]). Papy et al. verified that the Vandermonde
decomposition of $\mathcal{H}$ is

$\mathcal{H}=\mathcal{C}\times_{1}Z_{1}^{\top}\times_{2}Z_{2}^{\top}\dots\times_{m}Z_{m}^{\top},$

where $\mathcal{C}$ is a diagonal tensor with diagonal entries
$\\{c_{k}\\}_{k=1}^{K}$, and each $Z_{p}$ is a Vandermonde matrix

$Z_{p}^{\top}=\begin{bmatrix}1&z_{1}&z_{1}^{2}&\cdots&z_{1}^{I_{p}-1}\\\
1&z_{2}&z_{2}^{2}&\cdots&z_{2}^{I_{p}-1}\\\
\vdots&\vdots&\vdots&\vdots&\vdots\\\
1&z_{K}&z_{K}^{2}&\cdots&z_{K}^{I_{p}-1}\\\ \end{bmatrix}.$

So the target will be attained, if we obtain the Vandermonde decomposition of
this Hankel tensor $\mathcal{H}$. And the Vandermonde matrix can be estimated
by applying the total least square (TLS, cf. [10]) to the factor matrix in the
higher-order singular value decomposition (HOSVD, cf. [6, 12]) of the best
rank-$(K,K,\dots,K)$ approximation (cf. [7, 12]) of $\mathcal{H}$. Therefore,
computing the HOSVD of the best low rank approximation of a Hankel tensor is a
main part in exponential data fitting.

###  4.2 Fast HOOI for Hankel Tensors

Now we have a look at the best rank-$(R_{1},R_{2},\dots,R_{m})$ approximation
of a tensor. This concept was first introduced by De Lathauwer et al. in [7],
and they also proposed an effective algorithm called higher-order orthogonal
iterations (HOOI). There are other algorithms with faster convergence such as
[8] and [9] proposed, and one can refer to [12] for more details. However,
HOOI is still very popular, because it is so simple and still effective in
applications. Papy et al. chose HOOI in [18], so we aim at modifying the
general HOOI into a specific algorithm for Hankel tensors in order to get
higher speed for exponential data fitting. And we will focus on the square
Hankel tensor, since it is preferred in exponential data fitting.

The original HOOI algorithm for general tensors is as the following.

######  Algorithm 4.1.

HOOI for computing the best rank-$(R_{1},R_{2},\dots,R_{m})$ approximation of
tensor $\mathcal{A}\in\mathbb{C}^{I_{1}\times I_{2}\times\dots\times I_{m}}$.

Initialize $U_{p}\in\mathbb{C}^{I_{p}\times R_{p}}$ $(p=1:m)$ by HOSVD of
$\mathcal{A}$
Repeat
for $p=1:m$
$U_{p}\leftarrow R_{p}$ leading left singular vectors of
${\rm
Unfold}_{p}(\mathcal{A}\times_{1}\bar{U}_{1}\cdots\times_{p-1}\bar{U}_{p-1}\times_{p+1}\bar{U}_{p+1}\cdots\times_{m}\bar{U}_{m})$
end
Until convergence
$\mathcal{S}=\mathcal{A}\times_{1}\bar{U}_{1}\times_{2}\bar{U}_{2}\cdots\times_{m}\bar{U}_{m}$.

As employing this algorithm to Hankel tensors, we have to do some
modifications to involve the structures. First, the computation of HOSVD of a
Hankel tensor can be more efficient. The HOSVD is obtained by computing the
SVD of every unfolding of $\mathcal{A}$, Badeau and Boyer in [1] pointed out
that many columns in the unfolding of a structured tensor are redundant, so
that we can remove them for a smaller scale. Particularly, the mode-$p$
unfolding of a Hankel tensor $\mathcal{H}$ after removing the redundant
columns is a rectangular Hankel matrix of size
$I_{p}\times(d_{\mathcal{H}}-I_{p}+1)$ with the same generating vector with
$\mathcal{H}$

${\rm
Unfold}_{p}(\mathcal{H})\rightarrow\begin{bmatrix}h_{0}&h_{1}&\cdots&h_{d_{\mathcal{H}}-I_{p}-1}&h_{d_{\mathcal{H}}-I_{p}}\\\
h_{1}&\iddots&\iddots&\iddots&\vdots\\\
\vdots&\iddots&\iddots&\iddots&h_{d_{\mathcal{H}}-2}\\\
h_{I_{p}-1}&h_{I_{p}}&\cdots&h_{d_{\mathcal{H}}-2}&h_{d_{\mathcal{H}}-1}\end{bmatrix}.$

The SVD for Hankel matrices (called Symmetric SVD or Takagi Factorization) can
be computed more efficiently than for general matrices by using the algorithm
in [3, 26]. This is the first point that can be modified.

Next, there are plenty of Hankel tensor-matrix products in the above
algorithm, and they can be complemented by using our fast Hankel tensor-vector
products. For instance, the tensor-matrix product

$(\mathcal{H}\times_{2}\bar{U}_{2}\dots\times_{m}\bar{U}_{m})_{:,i_{2},\dots,i_{m}}=\mathcal{H}\times_{2}(\bar{U}_{2})_{:,i_{2}}\dots\times_{m}(\bar{U}_{m})_{:,i_{m}},$

and others are the same. So all the Hankel tensor-matrix products can be
computed by our fast algorithm.

Finally, if the Hankel tensor that we construct is further a square tensor,
then we can expect to obtain a _higher-order symmetric singular value
decomposition_ or called a _higher-order Takagi factorization_ (HOTF), that
is, $U_{1}=U_{2}=\dots=U_{m}$. Hence we do not need to update them separately
in each step. Therefore, the HOOI algorithm modified for square Hankel tensors
is as follows.

######  Algorithm 4.2.

HOOI for computing the best rank-$(R,R,\dots,R)$ approximation of square
Hankel tensor $\mathcal{H}\in\mathbb{C}^{I\times I\times\dots\times I}$.

$U\leftarrow R$ leading left singular vectors of ${\rm hankel}({\bf
h}_{0:I-1},{\bf h}_{I-1:nI-m})$
Repeat
$U\leftarrow R$ leading left singular vectors of ${\rm
Unfold}_{1}(\mathcal{H}\bar{U}^{m-1})$
Until convergence
$\mathcal{S}=\mathcal{H}\bar{U}^{m}$.

###  4.3 2-Dimensional Exponential Data Fitting

Many real-world signals are multi-dimensional, and the exponential model is
also useful for these problems, such that 2D curve fitting, 2D signal
recovery, etc. The multi-dimensional signals are certainly more complicated
than the one-dimensional signals. Thus one Hankel tensor is never enough for
representing the signals. So we need multi-level block Hankel tensors for
multi-dimensional signals. We take the 2D exponential data fitting for an
example to illustrate our block approach.

Assume that there is a 2D noiseless signal with $N_{1}\times N_{2}$ complex
samples $\\{x_{n_{1}n_{2}}\\}_{\begin{subarray}{l}n_{1}=0,1,\dots,N_{1}-1\\\
n_{2}=0,1,\dots,N_{2}-1\end{subarray}}$ which is modelled as a sum of $K$
exponential items

$x_{n_{1}n_{2}}=\sum_{k=1}^{K}a_{k}\exp({\bf
i}\varphi_{k})\exp\big{(}(-\alpha_{k}+{\bf i}\omega_{k})n_{1}\Delta
t_{1}+(-\beta_{k}+{\bf i}\nu_{k})n_{2}\Delta t_{2}\big{)},$

where the meanings of parameters are the same as those of 1D signals. Also,
the 2D signal can be rewritten into a compact form

$x_{n_{1}n_{2}}=\sum_{k=1}^{K}c_{k}z_{1,k}^{n_{1}}z_{2,k}^{n_{2}}.$

And our aim is still to estimate the poles $\\{z_{1,k}\\}$ and $\\{z_{2,k}\\}$
of the signal from the samples
$\\{x_{n_{1}n_{2}}\\}_{\begin{subarray}{l}n_{1}=0,1,\dots,N_{1}-1\\\
n_{2}=0,1,\dots,N_{2}-1\end{subarray}}$ that we receive.

Denote matrix $X=(x_{n_{1}n_{2}})_{N_{1}\times N_{2}}$. Then we construct a
BHHB tensor $\mathcal{H}$ with the generating matrix $X$ of a fixed order $m$,
size $(I_{1}J_{1})\times(I_{2}J_{2})\times\dots\times(I_{m}J_{m})$, and block
size $I_{1}\times I_{2}\times\dots\times I_{m}$. The order $m$ can be selected
arbitrarily and the size $I_{p}$ and $J_{p}$ of each dimension should be no
less than $K$. Then the BHHB tensor $\mathcal{H}$ has the level-2 Vandermonde
decomposition

$\mathcal{H}=\mathcal{C}\times_{1}(Z_{2,1}\oslash
Z_{1,1})^{\top}\times_{2}(Z_{2,2}\oslash
Z_{1,2})^{\top}\dots\times_{m}(Z_{2,m}\oslash Z_{1,m})^{\top},$

where $\mathcal{C}$ is a diagonal tensor with diagonal entries
$\\{c_{k}\\}_{k=1}^{K}$, each $Z_{1,p}$ or $Z_{2,p}$ is a Vandermonde matrix

$Z_{1,p}^{\top}=\begin{bmatrix}1&z_{1,1}&z_{1,1}^{2}&\cdots&z_{1,1}^{I_{p}-1}\\\
1&z_{1,2}&z_{1,2}^{2}&\cdots&z_{1,2}^{I_{p}-1}\\\
\vdots&\vdots&\vdots&\vdots&\vdots\\\
1&z_{1,K}&z_{1,K}^{2}&\cdots&z_{1,K}^{I_{p}-1}\\\ \end{bmatrix},\
Z_{2,p}^{\top}=\begin{bmatrix}1&z_{2,1}&z_{2,1}^{2}&\cdots&z_{2,1}^{J_{p}-1}\\\
1&z_{2,2}&z_{2,2}^{2}&\cdots&z_{2,2}^{J_{p}-1}\\\
\vdots&\vdots&\vdots&\vdots&\vdots\\\
1&z_{2,K}&z_{2,K}^{2}&\cdots&z_{2,K}^{J_{p}-1}\\\ \end{bmatrix},$

and the notation “$\oslash$” denotes the column-wise Kronecker product of two
matrices with the same column sizes, i.e.,

$(A\oslash B)(:,j)=A(:,j)\otimes B(:,j).$

So the target will be attained, if we obtain the level-2 Vandermonde
decomposition of this BHHB tensor $\mathcal{H}$.

We can use HOOI as well to compute the best rank-$(K,K,\dots,K)$ approximation
of the BHHB tensor $\mathcal{H}$

$\mathcal{H}=\mathcal{S}\times_{1}U_{1}^{\top}\times_{2}U_{2}^{\top}\dots\times_{m}U_{m}^{\top},$

where $\mathcal{S}\in\mathbb{C}^{K\times K\times\dots\times K}$ is the core
tensor and $U_{p}\in\mathbb{C}^{(I_{p}J_{p})\times K}$ has orthogonal columns.
Note that the tensor-vector products in HOOI should be computed by our fast
algorithm for BHHB tensor-vector products in Section 3. Then $U_{p}$ and
$Z_{2,p}\oslash Z_{1,p}$ have the common column space, which is equivalent to
that there is a nonsingular matrix $T$ such that

$Z_{2,p}\oslash Z_{1,p}=U_{p}T.$

Denote

$\begin{array}[]{l}A^{1\uparrow}=\big{[}A_{0:I-2,:}^{\top},A_{I,2I-2,:}^{\top},\dots,A_{(J-1)I:JI-2,:}^{\top}\big{]}^{\top},\\\
A^{1\downarrow}=\big{[}A_{1:I-1,:}^{\top},A_{I+1,2I-1,:}^{\top},\dots,A_{(J-1)I+1:JI-1,:}^{\top}\big{]}^{\top},\\\
A^{2\uparrow}=A_{0:(J-1)I-1,:},\\\ A^{2\downarrow}=A_{I:JI-1,:},\end{array}$

for matrix $A\in\mathbb{C}^{(IJ)\times K}$. Then it is easy to verify that

$\begin{split}(Z_{2,p}\oslash Z_{1,p})^{1\uparrow}D_{1}&=(Z_{2,p}\oslash
Z_{1,p})^{1\downarrow},\\\ (Z_{2,p}\oslash
Z_{1,p})^{2\uparrow}D_{2}&=(Z_{2,p}\oslash Z_{1,p})^{2\downarrow},\end{split}$

where $D_{1}$ is a diagonal matrix with diagonal entries
$\\{z_{1,k}\\}_{k=1}^{K}$ and $D_{2}$ is a diagonal matrix with diagonal
entries $\\{z_{2,k}\\}_{k=1}^{K}$. Then we have

$\begin{split}U_{p}^{1\uparrow}(TD_{1}T^{-1})&=U_{p}^{1\downarrow},\\\
U_{p}^{2\uparrow}(TD_{2}T^{-1})&=U_{p}^{2\downarrow}.\end{split}$

Therefore, if we obtain two matrices $W_{1}$ and $W_{2}$ satisfying that

$\begin{split}U_{p}^{1\uparrow}W_{1}&=U_{p}^{1\downarrow},\\\
U_{p}^{2\uparrow}W_{2}&=U_{p}^{2\downarrow},\end{split}$

then $W_{1}$ and $W_{2}$ share the same eigenvalues with $D_{1}$ and $D_{2}$,
respectively. Equivalently, the eigenvalues of $W_{1}$ and $W_{2}$ are exactly
the poles of the first and second dimension, respectively. Furthermore, we
also choose the total least square as in [18] for solving the above two
equations, since the noise on both sides should be taken into consideration.

######  Example 4.1.

We will show the effectiveness by a $2$-dimensional $2$-peak example partly
from [18],

$\begin{split}y_{n_{1}n_{2}}&=x_{n_{1}n_{2}}+e_{n_{1}n_{2}}\\\
&=\exp\big{(}(-0.01+2\pi{\bf
i}0.20)n_{1}\big{)}\cdot\exp\big{(}(-0.02+2\pi{\bf i}0.18)n_{2}\big{)}\\\
&\hskip 1.0pt+\exp\big{(}(-0.02+2\pi{\bf
i}0.22)n_{1}\big{)}\cdot\exp\big{(}(-0.01-2\pi{\bf
i}0.20)n_{2}\big{)}+e_{n_{1}n_{2}},\end{split}$

where the first dimension of this signal is the same as the example in [18]
and $e_{n_{1}n_{2}}$ is the noise. Since this signal has two peaks, the BHHB
tensor constructed following $\\{x_{n_{1}n_{2}}\\}$ is supposed to be rank-2.
Figure 3 illustrates the variation of the mode-1 singular values of a BHHB
tensor of size $6\times 6\times 6$ and block size $3\times 3\times 3$ with
different levels of noise. When the signal is noiseless, there are only two
nonzero mode-1 singular values, and others are about the machine epsilon, thus
is regarded as zero. When the noise is added on, the “zero singular values”
are at the same level with the noise, and those two largest singular values
can still be separated easily from the others.

Figure 3: The mode-1 singular values of the constructed BHHB tensor.

Then we compute the HOTF of the best rank-(2,2,2) approximation of a BHHB
tensor of size $30\times 30\times 30$ and block size $6\times 6\times 6$, and
apply the total least square twice to estimate the poles, as described
previously. Since we know the exact solution

$\begin{array}[]{ll}z_{11}=\exp(-0.01+2\pi{\bf
i}0.20),&z_{12}=\exp(-0.02+2\pi{\bf i}0.22),\\\ z_{21}=\exp(-0.02+2\pi{\bf
i}0.18),&z_{22}=\exp(-0.01-2\pi{\bf i}0.20),\end{array}$

the relative errors are shown in Figure 4. We conclude that the error can
attain the same level with the noise by employing our algorithm for 2D
exponential data fitting.

Figure 4: The relative errors with different levels of noise.

##  5 Conclusions

We propose a fast algorithm for Hankel tensor-vector products with
computational complexity $\mathcal{O}\big{(}(m+1)d_{\mathcal{H}}\log
d_{\mathcal{H}}\big{)}$, where $d_{\mathcal{H}}=n_{1}+n_{2}+\dots+n_{m}-m+1$.
This fast algorithm is derived by embedding the Hankel tensor into a larger
anti-circulant tensor, which can be diagonalized by the Fourier matrix. The
fast algorithm for block Hankel tensors with Hankel blocks is also described.
Furthermore, the fast Hankel and BHHB tensor-vector products are applied to
HOOI in 1D and 2D exponential data fitting, respectively. The numerical
experiments show the efficiency and effectiveness of our algorithms. Finally,
this fast scheme should be introduced into every algorithm which involves
Hankel tensor-vector products to improve its performance.

## Acknowledgements

Weiyang Ding would like to thank Prof. Sanzheng Qiao for the useful
discussions on fast algorithms for Hankel matrices. We also thank Professors
Lars Eldén and Michael K. Ng for their detailed comments.

## References

  * [1] R. Badeau and R. Boyer, Fast multilinear singular value decomposition for structured tensors, SIAM Journal on Matrix Analysis and Applications, 30 (2008), pp. 1008–1021.
  * [2] R. Boyer, L. De Lathauwer, K. Abed-Meraim, et al., Delayed exponential fitting by best tensor rank-$(R_{1},R_{2},R_{3})$ approximation, in IEEE Int. Conf. on Acoustics, Speech, and Signal Processing (ICASSP), 2005, pp. 19–23.
  * [3] K. Browne, S. Qiao, and Y. Wei, A Lanczos bidiagonalization algorithm for Hankel matrices, Linear Algebra and its Applications, 430 (2009), pp. 1531–1543.
  * [4] R. Chan and X.-Q. Jin, An Introduction to Iterative Toeplitz Solvers, SIAM, 2007.
  * [5] Z. Chen and L. Qi, Circulant tensors: Native eigenvalues and symmetrization, arXiv preprint arXiv:1312.2752, (2013).
  * [6] L. De Lathauwer, B. De Moor, and J. Vandewalle, A multilinear singular value decomposition, SIAM Journal on Matrix Analysis and Applications, 21 (2000), pp. 1253–1278.
  * [7] , On the best rank-1 and rank-$(R_{1},R_{2},...,R_{n})$ approximation of higher-order tensors, SIAM Journal on Matrix Analysis and Applications, 21 (2000), pp. 1324–1342.
  * [8] L. De Lathauwer and L. Hoegaerts, Rayleigh quotient iteration for the computation of the best rank-$(R_{1},R_{2},\dots,R_{N})$ approximation in multilinear algebra, tech. report, SCD-SISTA 04-003.
  * [9] L. Eldén and B. Savas, A Newton-Grassmann method for computing the best multilinear rank-$(r_{1},r_{2},r_{3})$ approximation of a tensor, SIAM Journal on Matrix Analysis and Applications, 31 (2009), pp. 248–271.
  * [10] G. H. Golub and C. F. Van Loan, Matrix Computations, The Johns Hopkins University Press, 4th ed., 2013.
  * [11] X.-Q. Jin, Developments and Applications of Block Toeplitz Iterative Solvers, Science Press, Beijing and Kluwer Academic Publishers, Dordrecht, 2002\.
  * [12] T. G. Kolda and B. W. Bader, Tensor decompositions and applications, SIAM Review, 51 (2009), pp. 455–500.
  * [13] L.-H. Lim, Singular values and eigenvalues of tensors: a variational approach, in Computational Advances in Multi-Sensor Adaptive Processing, 2005 1st IEEE International Workshop on, IEEE, 2005, pp. 129–132.
  * [14] F. T. Luk and S. Qiao, A fast singular value algorithm for Hankel matrices, in Contemporary Mathematics, vol. 323, Edited by Vadim Olshevsky, American Mathematical Society, Providence, RI; SIAM, Philadelphia, 2003, pp. 169–178.
  * [15] J.-G. Luque and J.-Y. Thibon, Hankel hyperdeterminants and Selberg integrals, Journal of Physics A: Mathematical and General, 36 (2003), pp. 5267–5292.
  * [16] M. K. Ng, Iterative Methods for Toeplitz Systems, Oxford University Press, 2004.
  * [17] V. Olshevsky, Structured Matrices in Mathematics, Computer Science, and Engineering. II, Contemporary Mathematics, vol. 281. American Mathematical Society, Providence, RI, 2001.
  * [18] J.-M. Papy, L. De Lathauwer, and S. Van Huffel, Exponential data fitting using multilinear algebra: The single-channel and multi-channel case, Numerical Linear Algebra with Applications, 12 (2005), pp. 809–826.
  * [19] , Exponential data fitting using multilinear algebra: The decimative case, Journal of Chemometrics, 23 (2009), pp. 341–351.
  * [20] V. Pereyra and G. Scherer, Exponential Data Fitting and its Applications, Bentham Science Publishers, 2010.
  * [21] L. Qi, Eigenvalues of a real supersymmetric tensor, Journal of Symbolic Computation, 40 (2005), pp. 1302–1324.
  * [22] , Eigenvalues and invariants of tensors, Journal of Mathematical Analysis and Applications, 325 (2007), pp. 1363–1377.
  * [23] , Hankel tensors: Associated Hankel matrices and Vandermonde decomposition, arXiv preprint arXiv:1310.5470, (2013).
  * [24] S. Qiao, G. Liu, and W. Xu, Block Lanczos tridiagonalization of complex symmetric matrices, in Advanced Signal Processing Algorithms, Architectures, and Implementations XV, Proceedings of the SPIE, vol. 5910, SPIE, 2005, pp. 285–295.
  * [25] Y. Song and L. Qi, Infinite and finite dimensional Hilbert tensors, arXiv preprint arXiv:1401.4966, (2014).
  * [26] W. Xu and S. Qiao, A fast symmetric SVD algorithm for square Hankel matrices, Linear Algebra and its Applications, 428 (2008), pp. 550–563.


# Neutrino parameters and the $N_{2}$-dominated scenario of
leptogenesis111Contribution to the Proceedings of the NuPhys2013 Conference,
19-20 December 2013, IOP, London.

M Re Fiorentina and SE Kinga,b a School of Physics and Astronomy, University
of Southampton, Southampton, SO17 1BJ, UK  b School of Physics and Astronomy,
Queen Mary University of London, London, E1 4NS, UK m.re-
<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

We briefly review the main aspects of leptogenesis, describing both the
unflavoured and the flavoured versions of the $N_{2}$-dominated scenario. A
study of the success rates of both classes of models has been carried out. We
comment on these results and discuss corrective effects to this simplest
scenario. Focusing on the flavoured case, we consider the conditions required
by strong thermal leptogenesis, where the final asymmetry is fully independent
of the initial conditions. Barring strong cancellations in the seesaw formula
and in the flavoured decay parameters, we show that strong thermal
leptogenesis favours a lightest neutrino mass $m_{1}\gtrsim 10$ meV for normal
ordering (NO) and $m_{1}\gtrsim 3$ meV for inverted ordering (IO). Finally, we
briefly comment on the power of absolute neutrino mass scale experiments to
either support or severely corner strong thermal leptogenesis.

##  1 Leptogenesis and the $N_{2}$-dominated scenario

Leptogenesis is a particularly attractive process for producing the baryon
asymmetry of the Universe, since it utilises the same mechanism that is able
to explain the observed neutrino masses and mixing, namely the seesaw
mechanism. Considering type-I seesaw, the Standard Model is extended by
including heavy right-handed (RH) Majorana neutrinos, $N_{i}$, that couple to
the lepton doublets via Yukawa interactions. Their $CP$\- and lepton-number-
violating decay produces a lepton asymmetry, which is partly converted into
the baryon sector by means of sphaleron processes. We assume in general the
presence of 3 RH neutrinos, thus giving rise to a model that depends on 18
parameters. These split into the 9 so-called low-energy neutrino parameters,
probed in experiments, and a further 9 describing the high-energy scale of the
RH neutrinos. We relate these two sets of parameters via the complex
orthogonal matrix $\Omega$ parameterisation [1]. By definition,
$N_{2}$-dominated [2] scenarios produce the observed asymmetry by the decay of
the next-to-lightest RH neutrino, $N_{2}$. The subsequent decay of the
lightest RH neutrino, $N_{1}$, produces a negligible contribution and,
furthermore, acts to wash the asymmetry out. These models are particularly
interesting because they naturally arise when conditions inspired by $SO(10)$
grand unification models are considered. Indeed, in this framework the
spectrum of the RH neutrinos is hierarchical with $M_{3}>10^{12}$ GeV,
$10^{9}\mbox{ GeV}<M_{2}<10^{12}\mbox{ GeV}$ and $M_{1}<10^{9}$ GeV, hence
leptogenesis is $N_{2}$-dominated. In the unflavoured case, flavour
interactions are neglected and RH neutrinos decay into a coherent
superposition of flavour eigenstates. The final value of the asymmetry is
exponentially suppressed by the unflavoured $N_{1}$ decay parameter,
$K_{1}\equiv\Gamma_{1}/H(T=M_{1})$, where $H$ is the Hubble parameter and
$\Gamma_{1}$ is the $N_{1}$ total decay width. Only when $K_{1}\ll 1$ is the
model able to produce the observed asymmetry. A numerical analysis shows that
this is realised in only $\sim 0.2\%$ of the parameter space.

This picture improves by taking account of flavour interactions. Considering
that $\tau$-interactions are efficient for $T\lesssim 10^{12}$ GeV, and muonic
ones for $T\lesssim 10^{9}$ GeV, $N_{2}$ decays in a two fully-flavoured
regime $10^{9}\mbox{ GeV}\lesssim T\lesssim 10^{12}\mbox{ GeV}$, while the
wash-out due to $N_{1}$ takes place in a three fully-flavoured regime
$T\lesssim 10^{9}$ GeV. Around the transition temperatures the behaviour
cannot be described by the usual Boltzmann equations and a density matrix
formalism must be adopted [3]. Avoiding these regions, the asymmetry produced
by $N_{2}$ is projected onto the flavour basis and the $N_{1}$ wash-out acts
separately on each flavour $\alpha$ ($\alpha=e,\,\mu,\,\tau$), exponentially
suppressing the asymmetry by the relative flavoured decay parameter,
$K_{1\alpha}\equiv(\Gamma_{i\alpha}+\bar{\Gamma}_{i\alpha})/H(T=M_{1})$. This
only requires one $K_{1\alpha}\ll 1$ in order to have a successful model. This
is found to happen in $\sim 30\%$ of the parameter space, showing that thanks
to flavour effects the $N_{2}$-dominated scenario can represent a viable model
of leptogenesis.

Some approximations are made: in particular, the presence of phantom terms and
flavour coupling is neglected [4]. The phantom terms correct the projection
onto the flavour basis of the asymmetry produced by $N_{2}$ in the coherent
superposition of $e$ and $\mu$ flavours. Each flavour asymmetry can be
regarded as a combination of a term proportional to the total lepton asymmetry
and a term related to the different flavour composition of the lepton quantum
states with respect to their $CP$ conjugate. This contribution is commonly
referred to as a phantom term. A further effect, flavour coupling, introduces
a correction when one considers that the Boltzmann equations are coupled due
to the asymmetries stored in the lepton doublets and in the Higgs bosons.
While the size of these effects has already been studied, we focus our
attention on the size of the parameter space region where these corrections
become significant. A numerical analysis reveals that for $\sim 20\%$ of the
parameter space the phantom terms and the flavour coupling modify the final
asymmetry by at least one order of magnitude.

##  2 Strong thermal leptogenesis

The final baryon asymmetry depends in general on the initial conditions. At
the high temperatures required by this minimal scenario, it is possible that
prior to leptogenesis some other mechanisms may generate a significant
asymmetry. Leptogenesis models able to wash out a pre-existing asymmetry
$N_{B-L}^{p,i}$, while producing the correct final amount, are said to satisfy
the strong thermal leptogenesis condition [5]. It has been shown [5] that this
case is realised only in a two-stage process where the RH neutrinos adhere to
the pattern discussed above, thus reinforcing the interest in
$N_{2}$-dominated models. In this scenario, $N_{2}$ decays in the two-
flavoured regime and efficiently washes out the $\tau$ component of
$N_{B-L}^{p,i}$. $N_{1}$, on the other hand, decays in the three-flavoured
regime, separately erasing the $e$ and $\mu$ components of $N_{B-L}^{p,i}$,
while the produced asymmetry survives in the $\tau$ flavour. This translates
into a set of conditions on the relevant decay parameters, namely
$K_{1e},K_{1\mu},K_{2\tau}\gg 1$ and $K_{1\tau}\lesssim 1$. When the low-
energy neutrino parameters are constrained within the current experimental
ranges, our analysis shows that, in the NO case, the strong thermal conditions
can be simultaneously satisfied, without fine-tuning, only if the absolute
neutrino mass scale $m_{1}$ is sufficiently large. In particular, for small
values of $m_{1}$ it becomes difficult to obtain a value of $K_{1e}$ large
enough to ensure the wash-out of a pre-existing asymmetry, while keeping
$K_{1\tau}$ small. It is then possible to place a precise analytical lower
bound on $m_{1}$ [6]. This lower bound depends on the magnitude of
$N_{B-L}^{p,i}$ and on the size of the entries of the orthogonal matrix
$\Omega$, i.e. on $\mbox{max}[\left|\Omega_{ij}\right|^{2}]\equiv M_{\Omega}$.
The larger $M_{\Omega}$, the more the seesaw mechanism relies, in order to
work, on fine-tuned cancellations, rather than on the natural interplay
between the electroweak scale and the higher, beyond the Standard Model, scale
of the RH neutrinos. In general, the lower bound shows a characteristic
dependence on the Dirac phase, $\delta$, and has its minimum for $\delta=0$.
For $M_{\Omega}\leq 2$ and a large pre-existing asymmetry,
$N_{B-L}^{p,i}=0.1$, at $\delta=0$ one gets $m_{1}\geq 0.7$ meV at $95\%$ C.L.
The lower bound is relaxed for lower values of $N_{B-L}^{p,i}$ and larger
values of $M_{\Omega}$. In particular, for $M_{\Omega}\gtrsim 4$ the lower
bound disappears. Our results are fully supported by a numerical analysis that
demonstrates the analytical lower bound being saturated at the expense of a
high level of fine-tuning in the seesaw formula and in the flavoured decay
parameters. Moreover, a statistical study of the $m_{1}$ distribution shows a
clear peak around $m_{1}\sim m_{\textup{\scriptsize atm}}\simeq 50$ meV and
rapid suppression for smaller values. In the aforementioned setup, more than
$99\%$ of points are found for $m_{1}\gtrsim 10$ meV (the value quoted in the
abstract).

In the IO case, tension exists between the smallness of $K_{1\tau}$ and a
$K_{1\mu}$ that must be large enough to wash-out the pre-existing asymmetry.
However, with respect to the NO case, $K_{1\mu}$ is easily large enough, even
for small values of $m_{1}$. More precisely, an analytical lower bound on
$m_{1}$ holds for IO only when $M_{\Omega}\lesssim 0.9$ and is therefore much
looser than for NO. The numerical analysis confirms again the analytical
results, but now $99\%$ of the points are found for $m_{1}\gtrsim 3$ meV (the
value quoted in the abstract).

It is important to stress again that these results rely on the need for a
$K_{1e}$ ($K_{1\mu}$), in NO (IO), large enough to efficiently wash-out the
pre-existing asymmetry, while keeping $K_{1\tau}$ small. This request becomes
nontrivial due to the experimental values of the neutrino mixing angles, in
particular to the suppression introduced by $\sin\theta_{13}$. It can be shown
that, with our usual requirements $N_{B-L}^{p,i}=0.1$ and $M_{\Omega}\leq 2$,
there is no lower bound on $m_{1}$ for $30^{\circ}\lesssim\theta_{13}\lesssim
50^{\circ}$. The reactor mixing angle happens to be large enough to allow both
itself and the Dirac phase to be probed, but not too large to trivially
satisfy strong thermal leptogenesis and remove any lower bound on $m_{1}$.

In the coming years, experimental results of the neutrino mass measurements,
especially from cosmology, may have the potential to test models of
leptogenesis. A measurement of the sum of the neutrino masses can put a
constraint on $m_{1}$ and, if the uncertainty is further reduced to $\sim 10$
meV [7], it will become possible to discriminate between vanishing and non-
vanishing $m_{1}$. Recent analyses that combine the data from Planck with the
latest growth of structure measurements from the Baryon Oscillation
Spectroscopic Survey (BOSS) have been able to give for the sum of the neutrino
masses a value $\sum_{i}m_{\nu_{i}}=\left(360\pm 100\right)$ meV [8], at
$3.4\sigma$ significance. Together with the measurements from neutrino
oscillation experiments this translates into $m_{1}=\left(127\pm 32\right)$
meV for NO and $m_{1}=\left(117\pm 32\right)$ meV for IO. These values are
quite high, but they serve as a strong hint at a non-vanishing absolute
neutrino mass scale, which goes in the right direction for strong thermal
leptogenesis. On the contrary, if new observation will measure $m_{1}<1$ meV,
a very little portion of the parameter space will let strong thermal
leptogenesis survive, at the expense of a high level of fine-tuning. Hence,
future cosmological observations have the power to either support or severely
corner strong thermal leptogenesis. However, this is highly dependent on the
ordering and therefore it is of the utmost importance that, in the next years,
neutrino oscillation experiments solve the ambiguity between NO and IO. The
former is in general a much more favourable case than IO for a significant
test, since it hints at larger deviations from the hierarchical limit
($m_{1}=0$).

##  3 Conclusions

Strong thermal leptogenesis introduces a theoretical request that is able to
constrain the parameter space and provide interesting results on $m_{1}$. In a
reasonable setup, favourable values of the absolute neutrino mass scale are
found for NO, $m_{1}\gtrsim 10$ meV and IO, $m_{1}\gtrsim 3$ meV. In NO a
further analytical lower bound is found for natural choices of $\Omega$
($M_{\Omega}\lesssim 4$), while in IO the analytical threshold is obtained
only for $M_{\Omega}\lesssim 0.9$. These constraints on $m_{1}$ allow the
future cosmological observations, together with the determination of the
neutrino mass ordering, to test strong thermal leptogenesis. We conclude by
highlighting that the link between cosmology and neutrino phenomenology can
prove to be once again extremely fertile and rich of ideas that may shed more
and more light on some of the major puzzles of modern physics.

## Acknowledgements

SK acknowledges financial support from the NExT/SEPnet Institute. MRF
acknowledges financial support from the STAG Institute. MRF is deeply grateful
to Pasquale Di Bari, Stefano Moretti and the School of Physics and Astronomy
of the University of Southampton. MRF wishes to thank Alexander J. Stuart and
is also profoundly indebted to Fiorenza Donato, Mariaelena Boglione and the
Physics Department of the University of Torino for the kind hospitality.

## References

## References

  * [1] Casas J A, Ibarra A, 2001 Nucl. Phys. B 618 171\.
  * [2] Di Bari P, 2005 Nucl. Phys. B 727 318-354.
  * [3] Blanchet S, Di Bari P, Jones D A, Marzola L, 2013 JCAP 1301 041
  * [4] Antusch S, Di Bari P, Jones D A, King S F, 2012 Nucl. Phys. B 856 180
  * [5] Bertuzzo E, Di Bari P, Marzola L, 2011 Nucl. Phys. B 849 521-548.
  * [6] Di Bari P, King S E, Re Fiorentin M, 2014 JCAP 1403 050
  * [7] Hamann J, Hannestad S, Wong Y, 2012 JCAP 1211 052
  * [8] Beutler F. et al. [BOSS Collaboration], arXiv:1403.4599 [astro-ph.CO].