#  Super-Mode OFDM With Index Modulation

Ali Tugberk Dogukan  and Ertugrul Basar  Manuscript received December 20,
2019; revised March 20, 2020; revised May 28, 2020; accepted July 17, 2020.
Date of publication XXXX XX, 2020; date of current version XXXX XX, 2020. This
work was supported by the Scientific and Technological Research Council of
Turkey (TÜBİTAK) under Grant 218E035. The associate editor coordinating the
review of this paper and approving it for publication was Prof. Zhaoyang
Zhang. The authors are with the Communications Research and Innovation
Laboratory (CoreLab), Department of Electrical and Electronics Engineering,
Koç University, Sariyer 34450, Istanbul, Turkey. Email<EMAIL_ADDRESS><EMAIL_ADDRESS>Color versions of one or more of the figures in this paper
are available online at http://ieeexplore.ieee.org. Digital Object Identifier
10.1109/TWC.2020.3010839

###### Abstract

Orthogonal frequency division multiplexing (OFDM) with index modulation (OFDM-
IM) appears as a promising multi-carrier waveform candidate for beyond 5G due
to its attractive advantages such as operational flexibility and ease of
implementation. However, OFDM-IM may not be a proper choice for 5G services
such as enhanced mobile broadband (eMBB) since achieving high data rates is
challenging because of its null subcarriers. One solution to enhance the
spectral efficiency of OFDM-IM is the employment of multiple distinguishable
constellations (modes) by also exploiting its null subcarriers for data
transmission. This paper proposes a novel IM technique called super-mode OFDM-
IM (SuM-OFDM-IM), where mode activation patterns (MAPs) and subcarrier
activation patterns (SAPs) are jointly selected and conventional data symbols
are repetition coded over multiple subcarriers to achieve a diversity gain.
For the proposed scheme, a low-complexity detector is designed, theoretical
analyses are performed and a bit error rate (BER) upper bound is derived. The
performance of the proposed system is also investigated through real-time
experiments using a software-defined radio (SDR) based prototype. We show that
SuM-OFDM-IM exhibits promising results in terms of spectral efficiency and
error performance; thus, appears as a potential candidate for 5G and beyond
communication systems.

###### Index Terms:

5G, index modulation (IM), subcarrier activation pattern (SAP), mode
activation pattern (MAP), diversity gain, log-likelihood ratio (LLR).

##  I Introduction

Index modulation (IM)[1], which employs the indices of transmit entities such
as antennas, time slots and radio frequency (RF) mirrors, has attracted
significant attention from the researchers due to its advantages over
conventional communication systems in terms of transceiver complexity and
spectral/energy efficiency. These attractive features of IM make it a
potential candidate for communication systems beyond 5G. IM has been employed
in numerous dimensions as mentioned above, and it can also be implemented in
the frequency domain by utilizing the indices of the available subcarriers of
multi-carrier systems[2, 3].

In [4], orthogonal frequency division multiplexing (OFDM) has been effectively
combined with IM by activating a number of subcarriers in a subblock while de-
activating the remaining ones. In OFDM-IM, information bits are conveyed by
not only $Q$-ary modulated symbols but also active subcarrier indices.
Compared to conventional OFDM, it is shown that OFDM-IM is capable of
providing a better error performance while consuming less power due to its
null subcarriers. The error performance of OFDM-IM can be improved by
employing a subcarrier-level block interleaver [5]. To further enhance the
error performance of OFDM-IM, coordinate interleaved OFDM-IM (CI-OFDM-IM),
which provides a diversity gain by transmitting the real and imaginary parts
of modulated symbols over different active subcarriers, is proposed in [6]. In
[7], OFDM with subcarrier number modulation (OFDM-SNM) is proposed where the
number of active subcarriers in a subblock is determined by the incoming bits.
A clever technique that decreases the modulation order while increasing the
number of bits transmitted by IM without any reduction in spectral efficiency,
is proposed [8]. An energy- and spectrum-efficient system that employs unused
subcarrier activation patterns (SAPs) in OFDM-IM by selecting a certain number
of adjacent subblocks jointly, is proposed in [9]. In recent years, OFDM-IM
has also been applied to emerging communication systems such as dual-hop
networks [10, 11], multiple-input multiple-output (MIMO) systems [12, 13, 14,
15], cognitive radio [16, 17] and optical communications [18].

Figure 1: Transmitter Structure of SuM-OFDM-IM.

Because of its null subcarriers, achieving very high data rates with OFDM-IM
is not possible and this may be critical for future enhanced mobile broadband
(eMBB) applications such as augmented/virtual reality (AR/VR) and $4$K/$8$K
video streaming. To overcome this critical drawback, multiple constellations
are utilized by activating not only the selected subcarriers but also the null
ones. Dual-mode aided OFDM (DM-OFDM), which also utilizes the null subcarriers
to transmit data symbols drawn from a secondary distinguishable constellation,
is proposed in [19]. In [20], a diversity gain is achieved by transmitting the
same data symbols through re-activation of the null subcarriers with a second
mode. DM-OFDM is generalized in [21] by changing the number of subcarriers
drawn from the same constellation. The performance of DM-OFDM and OFDM-IM has
been investigated in real-time by adopting software defined radios (SDRs) in
[22]. In DM-OFDM, only two different constellations are considered; however,
multiple-mode OFDM-IM (MM-OFDM-IM), which exploits all $n$ subcarriers in each
subblock with $n$ different distinguishable modes, is proposed to employ more
number of constellations [23]. MM-OFDM-IM uses the full permutations of these
constellations to apply IM. In [24], a generalized version of MM-OFDM-IM is
proposed where the size of the constellations can be altered although there is
no increase in the number of bits conveyed by IM. To obtain diversity gains,
coordinate interleaved MM-OFDM-IM (CI-MM-OFDM-IM) and linear constellation
precoded MM-OFDM-IM (LCP-MM-OFDM-IM) schemes have been presented [25].
Although these noticeable improvements, there is still an undeniable need for
multi-carrier waveforms that provide both reliable error performance and high
data rate for future eMBB services. For this reason, our motivation is to
exploit IM to further increase the spectral efficiency of OFDM-based waveforms
while providing a high reliability.

In this paper, we propose a novel OFDM-based transmission system called super-
mode OFDM-IM (SuM-OFDM-IM), which is capable of yielding both a diversity gain
and a high spectral efficiency. In SuM-OFDM-IM, according to the incoming
bits, mode activation patterns (MAPs) and SAPs are jointly determined;
therefore, additional information bits can be transmitted by IM. We also
consider determining MAPs and SAPs separately to provide an alternative
solution and call this scheme separated SuM-OFDM-IM (S-SuM-OFDM-IM). It is
well known that the bits transmitted by IM are more reliable than the ones
transmitted by conventional $Q$-ary PSK/QAM. In other words, conventional
PSK/QAM symbols restrict the diversity gain of the system. In our scheme, we
apply repetition coding to the modulated symbols drawn from the selected modes
over multiple subcarriers to attain a diversity gain. We propose a novel
encoding technique as well as a log-likelihood ratio (LLR)-based reduced
complexity maximum likelihood (ML) detection algorithm. Additionally, real-
time error performance of the proposed system is investigated through a
practical testbed using SDR units.

The rest of the paper is organized as follows. In Section II, we present the
system model of SuM-OFDM-IM. Section III deals with the selection of modes. In
Section IV, we derive an upper bound on the bit error rate (BER). Monte Carlo
simulation results and the complexity analysis are given in Section V. The
practical implementation of the proposed system is shown in Section VI and
finally, we conclude the paper in Section VII.

##  II System Model

In this section, we first introduce an encoding algorithm for SuM-OFDM-IM,
which jointly determines MAP and SAP in order to increase the number of bits
conveyed by IM. We also consider a clever symbol assignment strategy to
provide a diversity gain. Second, we propose a novel LLR-based reduced
complexity ML detector.

###  II-A Transmitter

SuM-OFDM-IM is an OFDM-based transmission system with $N$ subcarriers whose
transmitter structure is given in Fig. 1. A total of $m$ bits are conveyed for
each symbol transmission and split into groups with $p=m/g$ bits, where $g$ is
the number of subblocks. Since the same procedures are applied for each
subblock, for convenience, we consider the mapping of $p$ bits into the
$\alpha$th subblock $\mathbf{s}^{(\alpha)}$, where $\alpha\in\\{1,\dots,g\\}$.
These $p$ bits are split into two parts:

  1. 1.

$p_{1}=\lfloor{\log_{2}(\binom{M}{2}\binom{n}{n/2})}\rfloor$

  2. 2.

$p_{2}=\frac{n}{2}\log_{2}Q$

where $n$, $M$ and $Q$ are the size of the subblock,
$n=2^{r},r\in\\{2,3,\dots\\}$, the total number of modes, and constellation
size, respectively, and $\binom{\cdot}{\cdot}$ is the binomial coefficient.
First, $p_{1}$ bits jointly determine an MAP
$(\mathbf{v}^{(\alpha)}=[v_{1},v_{2}],v_{l}\in\\{1,\dots,M\\},l=1,2)$ and an
SAP
$(\mathbf{u}^{(\alpha)}=[u_{1},\dots,u_{n/2}],u_{m}\in\\{1,\dots,n\\},m=1,\dots,n/2)$,
where we discuss this joint selection in the sequel. Here, the selected MAP
includes two different modes to be able to employ IM $(v_{1}\neq v_{2})$ and
the strategy for mode selection is discussed in Section IV. We also activate
$n/2$ subcarriers per mode to maximize $p_{1}$. The remaining indices that are
not included in a selected SAP are represented by
$\mathbf{w}^{(\alpha)}=[u_{n/2+1},\dots,u_{n}]$. The subcarrier indices
$\mathbf{u}^{(\alpha)}$ and $\mathbf{w}^{(\alpha)}$ are exploited for the
symbols drawn from the first and second modes, respectively. Then, since each
symbol is repetition coded over a pair of subcarriers, $\mathbf{u}^{(\alpha)}$
and $\mathbf{w}^{(\alpha)}$ are split into $n/4$ groups
$(\mathbf{u}_{k}^{(\alpha)}=[u_{2k-1},u_{2k}],\mathbf{w}_{k}^{(\alpha)}=[u_{n/2+2k-1},u_{n/2+2k}],k=1,\dots,n/4$.
Second, $\frac{p_{2}}{2}=\frac{n}{4}\log_{2}Q$ bits determine the first $n/4$
modulated symbols $(\mathbf{x}_{1}^{(\alpha)}=[x_{1},\dots,x_{n/4}])$, which
are drawn from $\mathcal{M}_{v_{1}}$,
$\mathcal{M}_{i}=\\{\chi_{i1},\dots,\chi_{iQ}\\},i\in\\{1,\dots,M\\}$, where
$\chi_{iq}$ is the $q$th symbol of the $i$th mode and the set of all modes is
$\mathcal{M}=\\{\mathcal{M}_{1},\dots,\mathcal{M}_{M}\\}$. The average symbol
power of $\mathcal{M}$ is normalized to unity. Finally,
$\frac{p_{2}}{2}=\frac{n}{4}\log_{2}Q$ bits determine the last $n/4$ modulated
symbols $(\mathbf{x}_{2}^{(\alpha)}=[x_{n/4+1},\dots,x_{n/2}])$, which are
drawn from $\mathcal{M}_{v_{2}}$, and we obtain
$\mathbf{x}^{(\alpha)}=[\mathbf{x}^{(\alpha)}_{1}\mathbf{x}^{(\alpha)}_{2}]$.
The $k$th element of $\mathbf{x}_{1}^{(\alpha)}$ and
$\mathbf{x}_{2}^{(\alpha)}$ are repetition coded over the subcarriers of the
target subblock with the indices given by $\mathbf{u}_{k}^{(\alpha)}$ and
$\mathbf{w}_{k}^{(\alpha)}$, respectively. By performing the same steps for
all elements of $\mathbf{x}_{1}^{(\alpha)}$ and $\mathbf{x}_{2}^{(\alpha)}$,
$\mathbf{s}^{(\alpha)}$ is created. Consequently, a total of

$p=\Bigl{\lfloor}{\log_{2}\Bigl{(}\binom{M}{2}\binom{n}{n/2}\Bigl{)}}\Bigr{\rfloor}+\frac{n}{2}\log_{2}Q$
(1)

bits are transmitted per subblock.

The aim of this joint MAP and SAP selection is to further increase the number
of bits conveyed by IM. Our joint selection algorithm of MAP and SAP is
summarized below:

  1. 1.

Convert $p_{1}$ bits to a decimal number: $d$.

  2. 2.

Find the unique values $(a_{1},a_{2})$ satisfying the equation:

$d=a_{1}+\binom{M}{2}a_{2},$ (2)

where $a_{1}\in\\{0,\dots,\binom{M}{2}-1\\}$ and
$a_{2}\in\\{0,\dots,\binom{n}{n/2}-1\\}$ are the indices of the selected MAP
and SAP, respectively.

  3. 3.

Determine MAP $(\mathbf{v}^{(\alpha)})$ and SAP $(\mathbf{u}^{(\alpha)})$ by
applying $a_{1}$ and $a_{2}$ to the combinatorial algorithm [4], respectively.

TABLE I: Joint Selection of MAP and SAP $a_{1}$ |  $[v_{1},v_{2}]$ |  $a_{2}$ |  $[u_{1},u_{2}]$
---|---|---|---
$0$ |  $[1,2]$ |  $0$ |  $[1,2]$
$1$ |  $[1,3]$ |  $1$ |  $[1,3]$
$2$ |  $[2,3]$ |  $2$ |  $[2,3]$
$3$ |  $[1,4]$ |  $3$ |  $[1,4]$
$4$ |  $[2,4]$ |  $4$ |  $[2,4]$
$5$ |  $[3,4]$ |  $5$ |  $[3,4]$

Example: Assuming $M=4,n=4,Q=4$, the number of bits conveyed by IM is
$p_{1}=\lfloor{\log_{2}(\binom{4}{2}\binom{4}{4/2})}\rfloor=5$ and the number
of bits conveyed by conventional data symbols is
$p_{2}=\frac{4}{2}\log_{2}(4)=4$. For this case, a total of $\binom{4}{2}=6$
MAPs and $\binom{4}{2}=6$ SAPs are obtained as in Table I. There are $36$
possible pairs of $(a_{1},a_{2})$ and $32$ of these are employed to transmit
$\lfloor{\log_{2}(36)}\rfloor=5$ bits. Let $p_{1}$ and $p_{2}$ be
$\\{01001\\}$ and $\\{1110\\}$, respectively. The $(a_{1},a_{2})$ values
satisfying (2) are $(3,2)$. In the absence of joint selection, only four bits
can be transmitted per subblock instead of five by the indices of modes and
subcarriers for these given parameters. Therefore, an additional bit can be
conveyed by the joint selection in this specific case. Then, MAP is determined
as $\mathbf{v}=[1,4]$, i.e, $\mathcal{M}_{1}$ and $\mathcal{M}_{4}$ are
employed. SAP is determined as $\mathbf{u}=[1,3]$, so the remaining indices
are $\mathbf{w}=[2,4]$. For simplicity, we have removed the superscript
$(\alpha)$. The first two bits of $p_{2}$, which are $\\{11\\}$, select the
fourth symbol of the first selected mode $(\chi_{14})$ and the second two bits
of $p_{2}$, which are $\\{10\\}$, select the second symbol of the second
selected mode $(\chi_{42})$. Finally, $\chi_{14}$ and $\chi_{42}$ are placed
into the indices of subblock entries with $\mathbf{u}$ and $\mathbf{w}$,
respectively and the overall subblock is obtained as:

$\mathbf{s}=[\chi_{14}\ \chi_{42}\ \chi_{14}\ \chi_{42}]^{\mathrm{T}}.$ (3)

Having constructed the OFDM symbol subblock by subblock, a block-type
interleaver is executed as in [23] and the same procedures applied in OFDM are
performed. The inverse fast fourier transform (IFFT) algorithm is applied to
the OFDM block to obtain the time domain OFDM block as:

$\mathbf{s}_{T}=\mathrm{IFFT}\big{\\{}\mathbf{s}_{F}\big{\\}}=\begin{bmatrix}S(1)&S(2)&\cdots&S(N)\end{bmatrix}^{\mathrm{T}}.$
(4)

Then, a cyclic prefix (CP) of length $L$ samples
$\begin{bmatrix}S(N-L+1)&\cdots&S(N-1)S(N)\end{bmatrix}^{\mathrm{T}}$ is added
to the beginning of the OFDM block. After parallel to serial (P/S) and
digital/analog conversions are applied, the signal is transmitted over a
frequency-selective Rayleigh fading channel with the channel impulse response
(CIR) coefficients

$\mathbf{c}_{T}=\begin{bmatrix}c_{T}(1)&\cdots&c_{T}(v)\end{bmatrix}^{\mathrm{T}},$
(5)

where $c_{T}(\beta),\beta=1,\dots,v$ are circularly symmetric complex Gaussian
random variables with the $\mathcal{C}\mathcal{N}(0,\frac{1}{v})$
distribution. With the assumption that the channel remains constant during
transmission of an OFDM block and the CP length $L$ is larger than $v$, the
equivalent frequency domain input-output relationship of this OFDM scheme is
given by

$y_{F}(\beta)=s_{F}(\beta)c_{F}(\beta)+w_{F}(\beta),\hskip
5.69046pt\beta=1,\dots,N$ (6)

where $y_{F}(\beta)$, $c_{F}(\beta)$ and $w_{F}(\beta)$ are the received
signals, the channel fading coefficients and the noise samples in the
frequency domain, whose vector forms are given as $\mathbf{y}_{F}$,
$\mathbf{c}_{F}$ and $\mathbf{w}_{F}$, respectively. The distributions of
$c_{F}(\beta)$ and $w_{F}(\beta)$ are $\mathcal{C}\mathcal{N}(0,1)$ and
$\mathcal{C}\mathcal{N}(0,N_{0})$, respectively, where $N_{0}$ is the noise
variance in the frequency domain, which is equal to the noise variance in the
time domain. The signal-to-noise ratio (SNR) is defined as $\rho=E_{b}/N_{0}$,
where $E_{b}=(N+L)/m$ is the average transmitted energy per bit. For
simplicity, the effect of CP on the spectral efficiency is ignored. Then, the
spectral efficiency of the proposed scheme is obtained as

$\eta=m/N\hskip 5.69046pt\text{[bits/s/Hz]}.$ (7)

####  II-A1 Separate selection of MAPs and SAPs

Instead of joint selection, it is possible to determine MAPs and SAPs,
separately. This scheme is called S-SuM-OFDM-IM. The only difference from SuM-
OFDM-IM is the selection procedure of MAPs and SAPs. In this case, a total
number of

$p=\Big{\lfloor}{\log_{2}\Big{(}\binom{M}{2}\Big{)}}\Big{\rfloor}+\Big{\lfloor}{\log_{2}\Big{(}\binom{n}{n/2}\Big{)}}\Big{\rfloor}+\frac{n}{2}\log_{2}Q$
(8)

bits are transmitted per subblock. As seen from Fig. 2, the joint selection
provides higher spectral efficiency than the separate one for varying $Q$
values. However, in Section IV, we provide an interesting trade-off between
error performance and spectral efficiency.

Figure 2: Spectral efficiency comparison of joint and separate selections of
MAPs and SAPs for varying values of $Q$ when $n=4$ and $M=4$.

###  II-B Receiver

Algorithm 1 LLR-Based Reduced Complexity ML Detector $(n\in\\{4,8\\})$

1:$\mathbf{y}_{F}^{(\alpha)}$, $\mathbf{c}_{F}^{(\alpha)}$, $N_{0}$,
$\mathcal{M}$, $M$, $Q$, $p_{1}$, $n$

2:$\hat{a_{1}},\hat{a_{2}},\hat{\mathbf{x}}$

3:$\mathbf{y}:=\mathbf{y}_{F}^{(\alpha)},\mathbf{c}:=\mathbf{c}_{F}^{(\alpha)}$

4:for $i=1$ to $\binom{n}{n/2}$ do

5: $\mathcal{I}^{(i)}:=$combinatorics_n2c($i-1$,$n/2$)

6: $j:=1$

7: for $m=1$ to $M$ do

8: for $h=1$ to $n/4$ do

9: for $q=1$ to $Q$ do

10: for $z=1$ to $2$ do

11: $\hat{\mathcal{I}}:=\mathcal{I}^{(i,z+2(h-1))}+1$

12:
$\delta_{1}^{(q,z+2(h-1))}:=-|\mathbf{y}^{(\hat{\mathcal{I}})}-\mathbf{c}^{(\hat{\mathcal{I}})}\mathcal{M}^{(m,q)}|^{2}/N_{0}$

13: $\delta_{2}^{(q,h)}:=\delta_{1}^{(q,2h-1)}+\delta_{1}^{(q,2h)}$

14: end for

15: end for

16: end for

17: for $g_{1}=1$ to $n/4$ do

18: for $g_{2}=1$ to $g_{1}$ do

19: if $n=8$ $\&\&$ $g_{1}\neq g_{2}$ then

20: for $q_{1}=1$ to $Q$ do

21: for $q_{2}=1$ to $Q$ do

22: $\gamma^{(j,i)}:=\delta_{2}^{(q_{1},g_{1})}+\delta_{2}^{(q_{2},g_{2})}$

23: $j:=j+1$

24: end for

25: end for

26: else if $n=4$ then

27: for $q_{1}=1$ to $Q$ do

28: for $q_{2}=1$ to $Q$ do

29: $\gamma^{(j,i)}:=\delta_{2}^{(q_{1})}+\delta_{2}^{(q_{2})}$

30: $j:=j+1$

31: end for

32: end for

33: end if

34: end for

35: end for

36: end for

37:end for

38:for $i=1$ to $\binom{n}{n/2}$ do

39: $j:=1$

40: for $m=1$ to $MQ^{2}$ with increments $Q^{2}$ do

41: $\Delta:=$MAX($\gamma^{(m,i)},\gamma^{(m+1,i)}$)

42: for $q=2$ to $Q^{2}-1$ do

43: $\Delta:=$MAX($\Delta,\gamma^{(m+q,i)}$)

44: end for

45: $\Gamma^{(j,i)}:=\Delta$

46: for $z=1$ to $n/2$ do

47: $\hat{\mathcal{I}}:=\mathcal{I}^{(i,z)}$

48: $\epsilon^{(z)}:=-|\mathbf{y}^{(\hat{\mathcal{I}})}|^{2}/N_{0}$

49: end for

50: $\Gamma^{(j,i)}:=\Gamma^{(j,i)}+\sum_{z=1}^{n/2}\epsilon^{(z)}$

51: $j:=j+1$

52: end for

53:end for

54:$\Gamma:=\Gamma^{\mathrm{T}}$

53:for $v=1$ to $\binom{M}{2}$ do

54: $\tilde{\mathcal{I}}:=$combinatorics_n2c($v-1,2$) $+1$

55: for $r=1$ to $\binom{n}{n/2}$ do

56:
$\mathbf{\Lambda}^{(r,v)}:=\Gamma^{(r,\tilde{\mathcal{I}}^{(1)})}+\Gamma^{(\binom{n}{n/2}-r+1,\tilde{\mathcal{I}}^{(2)})}$

57: end for

58:end for

59:for $a_{1}=0$ to $\binom{M}{2}-1$ do

60: for $a_{2}=0$ to $\binom{n}{n/2}-1$ do

61: if ($a_{1}+\binom{M}{2}a_{2}\geq 2^{p_{1}}$) then

62: $\mathbf{\Lambda}^{(a_{2}+1,a_{1}+1)}:=-\infty$

63: end if

64: end for

65:end for

66:$(\hat{a_{2}},\hat{a_{1}})=\underset{a_{1},a_{2}}{\arg\>\max}\hskip
2.84544pt{\mathbf{\Lambda}^{(a_{2},a_{1})}}$

67:$\hat{\mathbf{v}}:=$$\mathrm{combinatorics\\_n2c}$($\hat{a_{1}}-1$,$2$)$+1$

68:$\hat{\mathbf{u}}:=$$\mathrm{combinatorics\\_n2c}$($\hat{a_{2}}-1$,$n/2$)$+1$

69:$\hat{\mathbf{w}}:=$$\mathrm{combinatorics\\_n2c}$($\binom{n}{n/2}-\hat{a_{2}}-1$,$n/2$)$+1$

70:if $n=8$ then

71:
$\hat{x}_{1}:=\underset{q}{\arg\>\min}{\left\|[\mathbf{y}^{(\hat{\mathbf{u}}^{(1)})}\mathbf{y}^{(\hat{\mathbf{u}}^{(2)})}]-[\mathbf{c}^{(\hat{\mathbf{u}}^{(1)})}\mathbf{c}^{(\hat{\mathbf{u}}^{(2)})}]\mathcal{M}^{(\hat{\mathbf{v}}^{(1)},q)}\right\|}^{2}$

72:
$\hat{x}_{2}:=\underset{q}{\arg\>\min}{\left\|[\mathbf{y}^{(\hat{\mathbf{u}}^{(3)})}\mathbf{y}^{(\hat{\mathbf{u}}^{(4)})}]-[\mathbf{c}^{(\hat{\mathbf{u}}^{(3)})}\mathbf{c}^{(\hat{\mathbf{u}}^{(4)})}]\mathcal{M}^{(\hat{\mathbf{v}}^{(1)},q)}\right\|}^{2}$

73:
$\hat{x}_{3}:=\underset{q}{\arg\>\min}{\left\|[\mathbf{y}^{(\hat{\mathbf{w}}^{(1)})}\mathbf{y}^{(\hat{\mathbf{w}}^{(2)})}]-[\mathbf{c}^{(\hat{\mathbf{w}}^{(1)})}\mathbf{c}^{(\hat{\mathbf{w}}^{(2)})}]\mathcal{M}^{(\hat{\mathbf{v}}^{(2)},q)}\right\|}^{2}$

74:
$\hat{x}_{4}:=\underset{q}{\arg\>\min}{\left\|[\mathbf{y}^{(\hat{\mathbf{w}}^{(3)})}\mathbf{y}^{(\hat{\mathbf{w}}^{(4)})}]-[\mathbf{c}^{(\hat{\mathbf{w}}^{(3)})}\mathbf{c}^{(\hat{\mathbf{w}}^{(4)})}]\mathcal{M}^{(\hat{\mathbf{v}}^{(2)},q)}\right\|}^{2}$

75: $\hat{\mathbf{x}}=[\hat{x}_{1}\hskip 2.84544pt\hat{x}_{2}\hskip
2.84544pt\hat{x}_{3}\hskip 2.84544pt\hat{x}_{4}]$

76:else if $n=4$ then

77:
$\hat{x}_{1}:=\underset{q}{\arg\>\min}{\left\|[\mathbf{y}^{(\hat{\mathbf{u}}^{(1)})}\mathbf{y}^{(\hat{\mathbf{u}}^{(2)})}]-[\mathbf{c}^{(\hat{\mathbf{u}}^{(1)})}\mathbf{c}^{(\hat{\mathbf{u}}^{(2)})}]\mathcal{M}^{(\hat{\mathbf{v}}^{(1)},q)}\right\|}^{2}$

78:
$\hat{x}_{2}:=\underset{q}{\arg\>\min}{\left\|[\mathbf{y}^{(\hat{\mathbf{w}}^{(1)})}\mathbf{y}^{(\hat{\mathbf{w}}^{(2)})}]-[\mathbf{c}^{(\hat{\mathbf{w}}^{(1)})}\mathbf{c}^{(\hat{\mathbf{w}}^{(2)})}]\mathcal{M}^{(\hat{\mathbf{v}}^{(2)},q)}\right\|}^{2}$

79: $\hat{\mathbf{x}}=[\hat{x}_{1}\hskip 2.84544pt\hat{x}_{2}]$

80:end if
return $\hat{a_{1}},\hat{a_{2}},\hat{\mathbf{x}}$

At the receiver side, after removing the CP, de-interleaving and performing
FFT, ML detection can be considered to decode the $\alpha$th subblock. The
subblock set, which includes all possible subblock realizations, is defined as

$\mathcal{S}=\\{\mathbf{s}_{1},\mathbf{s}_{2},\dots\mathbf{s}_{2^{p}}\\}.$ (9)

Then, ML detection rule for SuM-OFDM-IM is given as

$\hat{\mathbf{s}}=\underset{\mathbf{s}\in\mathcal{S}}{\arg\>\min}{\left\|\mathbf{y}_{F}^{(\alpha)}-\mathbf{c}_{F}^{(\alpha)}\odot\mathbf{s}\right\|}^{2},$
(10)

where $\mathbf{y}_{F}^{(\alpha)}$ and $\mathbf{c}_{F}^{(\alpha)}$ and are the
vectors of the received signals and the channel fading coefficients
corresponding to $\alpha$th subblock, respectively, and $\odot$ is the
element-wise product.

The ML detection complexity per subblock is proportional to $2^{p}$,
therefore, it is not practical and efficient for large values of $n$, $M$ and
$Q$. Thus, we propose a novel LLR-based reduced complexity ML detection
algorithm for $n=4$ and $n=8$ as seen in Algorithm 1, which significantly
reduces the detection complexity as discussed in Section V. Since we perform
repetition coding for data symbols, we present a modified LLR detector. This
detector achieves the same performance as the ML detector since it searches
for all possible MAP and SAP combinations.

The functions considered in Algorithm 1 are
$\mathrm{MAX}(a,b)=\mathrm{max}(a,b)+\ln(1+e^{-|a-b|})$ and the combinatorial
algorithm $\mathcal{I}=\mathrm{combinatorics\\_n2c}(c,d)$ that takes the
integer $c$ and the number of indices $d$ as inputs and outputs the set of
indices $\mathcal{I}$ [4].

The major steps of Algorithm 1 are explained below:

  1. 1.

The inputs of the algorithm are $\alpha$th received subblock
$(\mathbf{y}_{F}^{(\alpha)})$, the channel fading coefficients corresponding
to $\alpha$th subblock $(\mathbf{c}_{F}^{(\alpha)})$, noise power $(N_{0})$,
the set of all modes $(\mathcal{M})$, the number of modes $(M)$, modulation
order $(Q)$, the number of bits conveyed by IM $(p_{1})$ and the subblock size
$(n)$.

  2. 2.

To determine the most likely MAP and SAP combination, LLR values are obtained
between the lines (2-58). Specifically, between the lines (2-35), we provide a
clever LLR calculation technique due to the unique symbol assignment strategy
of our scheme. An iterative approach is applied between the lines (39-42) as
in [19].

  3. 3.

By assigning $-\xi$, $(\xi>>1)$ to LLR values corresponding to illegal MAP and
SAP combinations between the lines (59-65), we avoid them. Thus, these illegal
values can not be selected.

  4. 4.

Based on LLR values, MAP and SAP are jointly detected between the lines
(66-69). For example, assuming $M=4,n=4$, we obtain

$\mathbf{\Lambda}=\begin{bmatrix}\Lambda_{11}&\Lambda_{12}&\Lambda_{13}&\Lambda_{14}&\Lambda_{15}&\Lambda_{16}\\\
\Lambda_{21}&\Lambda_{22}&\Lambda_{23}&\Lambda_{24}&\Lambda_{25}&\Lambda_{26}\\\
\Lambda_{31}&\Lambda_{32}&\Lambda_{33}&\Lambda_{34}&\Lambda_{35}&\Lambda_{36}\\\
\Lambda_{41}&\Lambda_{42}&\Lambda_{43}&\Lambda_{44}&\Lambda_{45}&\Lambda_{46}\\\
\Lambda_{51}&\Lambda_{52}&\Lambda_{53}&\Lambda_{54}&\Lambda_{55}&\Lambda_{56}\\\
\Lambda_{61}&\Lambda_{62}&-\xi&-\xi&-\xi&-\xi\\\ \end{bmatrix},$ (11)

where $\Lambda_{ij}$ is the LLR value in $i$th row and $j$th column of the LLR
matrix $\mathbf{\Lambda}$, which is formed between the lines (53-58). Columns
and rows of $\mathbf{\Lambda}$ represents the possible MAPs and SAPs,
respectively. As seen from (11), $\mathbf{\Lambda}$ includes LLR values for
all possible MAP and SAP combinations; therefore, this detector provides the
same error performance as the ML detector.

  5. 5.

After detecting the active MAP and SAP, data symbols are decoded between the
lines (70-80).

  6. 6.

Finally, algorithm outputs the decoded data symbols $(\hat{\mathbf{x}})$ and
decimal equivalent of IM bits $(\hat{d}=\hat{a}_{1}+\binom{M}{2}\hat{a}_{2})$.

Remark (Generalization): To generalize the LLR-based detector for arbitrary
$n$, only the following changes should be performed in Algorithm 1:

  * •

The condition in line 17 is changed as $n\neq 4$ rather than $n=8$.

  * •

The lines between (70-80) are modified. For any values of $n$, a total of
$n/2$ metric calculations are performed. If $\kappa\leq n/4$, the parameters
$\hat{u}^{(2\kappa-1)}$, $\hat{u}^{(2\kappa)}$, and $\hat{v}^{(1)}$ are
exploited to detect $\hat{x}_{\kappa}$, $\kappa=1,\cdots,n/2$. If
$\kappa>n/4$, the parameters $\hat{w}^{(2(\kappa-n/4)-1)}$,
$\hat{w}^{(2(\kappa-n/4))}$, and $\hat{v}^{(2)}$ are exploited to detect
$\hat{x}_{\kappa}$.

##  III Mode Selection for SuM-OFDM-IM

In this section, we discuss our strategy for mode selection. According to
[23], the optimal modes must maximize the minimum intra-mode distance (MIAD)

$\displaystyle\epsilon_{1}$
$\displaystyle=\underset{\iota,\kappa\in\\{1,\dots,2^{p}\\}}{\min}{\left\|\mathbf{S}_{\iota}-\mathbf{S}_{\kappa}\right\|}^{2}_{F},$
s.t. $\displaystyle\hskip
5.69046ptE\\{\|\mathbf{S}\|^{2}_{F}\\}=n\>\>\text{and}\>\>\mathrm{rank}(\mathbf{S}_{\iota}-\mathbf{S}_{\kappa})=1,$
(12)

where $\mathbf{S}=\mathrm{diag}(\mathbf{s})$, $\mathrm{diag}(\mathbf{a})$
represents a diagonal matrix whose diagonal elements are $\mathbf{a}$, and the
minimum inter-mode distance (MIRD)

Figure 3: Partition of $16$-QAM constellation for $M=4$ and $Q=4$.

$\displaystyle\epsilon_{2}$
$\displaystyle=\underset{\iota,\kappa\in\\{1,\dots,2^{p}\\}}{\min}{\left\|\mathbf{S}_{\iota}-\mathbf{S}_{\kappa}\right\|}^{2}_{F},$
s.t. $\displaystyle\hskip
5.69046ptE\\{\|\mathbf{S}\|^{2}_{F}\\}=n\>\>\text{and}\>\>\mathrm{rank}(\mathbf{S}_{\iota}-\mathbf{S}_{\kappa})=2,$
(13)

where $\mathbf{S}_{\iota}$ and $\mathbf{S}_{\kappa}$ are two different
realizations of $\mathbf{S}$. For simplicity, we have removed superscript
$(\alpha)$ in $\mathbf{S}^{(\alpha)}$, since all subblocks are identical. The
modes can be created by partitioning PSK/QAM constellations as in [23] and it
is proven that QAM constraint performs better than PSK constraint for $MQ>4$.
Therefore, we focus on only QAM constellation for our system model:

  * •

QAM constraint to maximize MIRD:

$d_{\mathrm{inter}}(M,Q)=\begin{cases}2\sqrt{\frac{6}{5QM-4}},&\text{for
rectangular $MQ-$QAM}\\\ \sqrt{\frac{6}{QM-1}},&\text{for square
$MQ-$QAM}\end{cases}$ (14)

  * •

QAM constraint to maximize MIAD:

$\displaystyle d_{\mathrm{intra}}(M,Q)=$
$\displaystyle\begin{cases}d_{\mathrm{inter}}(M,Q)\frac{\sqrt{5M}}{2},&\text{for
rectangular $MQ-$QAM ($Q=2$)}\\\
d_{\mathrm{inter}}(M,Q)\sqrt{M},&\text{otherwise}\end{cases}$ (15)

  * •

QAM constraint to maximize MIAD for $MQ>>1$:

$\displaystyle d_{\mathrm{intra}}(M,Q)\approx d_{\mathrm{intra}}(Q)=$
$\displaystyle\begin{cases}2\sqrt{\frac{6}{5Q}},&\text{for rectangular
$MQ-$QAM ($Q\neq 2$)}\\\ \sqrt{\frac{6}{Q}},&\text{otherwise}\end{cases}$ (16)

with these constraints, modes can be obtained by partitioning QAM
constellations for given parameters $(M,Q)$. For example, for $M=4,Q=4$, the
optimal modes for $16$-QAM can be selected as in Fig. 3, where
$d_{\mathrm{inter}}=0.6325$ and $d_{\mathrm{intra}}=1.2649$ from (14) and
(15), respectively.

##  IV Performance Analysis

In this section, we derive an upper bound on the BER of the SuM-OFDM-IM system
employing ML detection. From (10), the conditional pairwise error probability
(CPEP), which is defined as the probability of transmitting $\mathbf{S}$ and
erroneously detecting $\hat{\mathbf{S}}$ conditioned on $\mathbf{c}_{F}$, is
given as:

$\text{Pr}(\mathbf{S}\rightarrow\hat{\mathbf{S}}|\mathbf{c}_{F})=\mathcal{Q}\bigg{(}\sqrt{\frac{\rho}{2}\Big{\|}(\mathbf{S}-\hat{\mathbf{S}})\mathbf{c}_{F}\Big{\|}^{2}}\bigg{)}.$
(17)

Then, by approximating
$\mathcal{Q}(x)\approx\frac{e^{-x^{2}/2}}{12}+\frac{e^{-2x^{2}/3}}{4}$, the
unconditioned PEP (UPEP) can be obtained as in [4]:

$\displaystyle\text{Pr}\big{(}\mathbf{S}\rightarrow\hat{\mathbf{S}}\big{)}$
$\displaystyle=E_{\mathbf{c}_{F}}\\{\text{Pr}(\mathbf{S}\rightarrow\hat{\mathbf{S}}|\mathbf{c}_{F})\\}$
$\displaystyle\approx\frac{1/12}{\text{det}(\mathbf{I}_{n}+\rho_{1}\mathbf{A})}+\frac{1/4}{\text{det}(\mathbf{I}_{n}+\rho_{2}\mathbf{A})},$
(18)

where $\rho_{1}=1/(4N_{0}),\rho_{2}=1/(3N_{0})$,
$\mathbf{A}=(\mathbf{S}-\hat{\mathbf{S}})^{\mathrm{H}}(\mathbf{S}-\hat{\mathbf{S}})$
and $(.)^{\mathrm{H}}$ denotes Hermitian transposition. Finally, according to
the union bounding technique, the BER of SuM-OFDM-IM can be upper bounded by

$P_{e}\leq\frac{1}{p2^{p}}\sum_{\mathbf{S}}\sum_{\hat{\mathbf{S}}}\text{Pr}\big{(}\mathbf{S}\rightarrow\hat{\mathbf{S}}\big{)}e\big{(}\mathbf{S},\hat{\mathbf{S}}\big{)},$
(19)

where Pr$(\mathbf{S}\rightarrow\hat{\mathbf{S}})$ is given by (18) and
$e(\mathbf{S},\hat{\mathbf{S}})$ represents the number of bit errors for the
corresponding pairwise error event.

We have
$\text{det}(\mathbf{A})=\prod_{\zeta=1}^{r}\lambda_{\zeta}(\mathbf{A})$ where
$\lambda_{\zeta}(\mathbf{A})$ is the $\zeta$th eigenvalue of $\mathbf{A}$ and
$r=\mathrm{rank}(\mathbf{A})$ determines the diversity order of the system,
since $P_{e}$ decays with $r$. At high SNR values, by ignoring
$\mathbf{I}_{n}$ terms in (18), (19) can be simplified as follows:

$P_{e}\leq\frac{\rho_{1}^{-r}+3\rho_{2}^{-r}}{12p2^{p}}\sum_{\mathbf{S}}\sum_{\hat{\mathbf{S}}}\Big{(}\prod_{\zeta=1}^{r}\lambda_{\zeta}(\mathbf{A})\Big{)}^{-1}e\big{(}\mathbf{S},\hat{\mathbf{S}}\big{)}.$
(20)

We investigate two cases: $(i)$ erroneous detection of index bits and $(ii)$
erroneous detection of a single or multiple data symbols while the index bits
are decoded correctly. For case $(i)$, due to our novel MAP and SAP selection
strategy, we always obtain $r\geq 2$. For example, from Table I, assume that
$\mathbf{v}=[1,2]$ is selected and $\hat{\mathbf{v}}=[1,3]$ is decoded
erroneously, which is the worst case error event due to overlapping index
values of $\mathbf{v}$ and $\hat{\mathbf{v}}$. Even in this error event, $r$
is obtained as $2$ as in plain OFDM-IM systems [4]. If $\mathbf{v}$ and
$\hat{\mathbf{v}}$ have non-overlapping indices, $r$ might be even equal to
$3$ or $4$. To explain case $(ii)$, we illustrate the following example.
Assume that $n=4$, $M=4$, $Q=4$ and
$\mathbf{S}=\mathrm{diag}([\chi_{11}\>\chi_{11}\>\chi_{31}\>\chi_{31}])$,
which is obtained as in (3), is transmitted and
$\hat{\mathbf{S}}=\mathrm{diag}([\chi_{11}\>\chi_{11}\>\chi_{32}\>\chi_{32}])$
is erroneously detected, i.e, MAP and SAP are detected correctly; however, the
second symbol is decoded erroneously $(\chi_{31}\rightarrow\chi_{32})$. Then,
we obtain:

$\mathbf{A}=\begin{bmatrix}0&0&0&0\\\ 0&0&0&0\\\ 0&0&1.6&0\\\
0&0&0&1.6\end{bmatrix},$ (21)

where $r=\mathrm{rank}(\mathbf{A})=2$. As seen from this example, due to
clever symbol assignment strategy of our scheme, even the erroneous detection
of a single data symbol ensures a diversity order of $2$. It can be easily
shown that for all error events of this type, we again obtain at least $r=2$.

TABLE II: Percentage of Error Events $(M,Q)$ |  $r=2$ |  $r=3$ |  $r=4$
---|---|---|---
$(4,4)$ |  $4.79\%$ |  $15.07\%$ |  $80.14\%$
$(8,2)$ |  $5.10\%$ |  $14.95\%$ |  $79.95\%$
$(4,16)$ |  $1.14\%$ |  $4.02\%$ |  $94.84\%$
$(16,4)$ |  $1.27\%$ |  $3.92\%$ |  $94.81\%$

To provide further insights on the performance of SuM-OFDM-IM, we obtain the
percentage of error events with $r=2,3$ and $4$ in Table II for varying
$(M,Q)$ by considering all possible error events for $n=4$. It should be noted
that the percentage of error events has a remarkable affect on the overall
error performance of a system from (19). As it can be deduced from Table II,
since $r=4$ is obtained much more frequently than other possible rank values,
our proposed system is capable of providing a superior error performance. It
can be seen from Table II that as the overall constellation size $(MQ)$
increases, the percentage of $r=4$ increases due to the increase in the number
of data symbols. Because, when MAP and SAP are detected correctly and two data
symbols are decoded erroneously, we obtain $r=4$. Finally, in light of the
above discussion, we obtain

$\underset{\mathbf{S},\hat{\mathbf{S}}}{\min}\ \mathrm{rank}(\mathbf{A})=2,$
(22)

which proves that our system provides a diversity order of $2$.

##  V Simulation Results

Figure 4: Theoretical BER upper bound of SuM-OFDM-IM with simulation results.

In this section, we perform Monte Carlo simulations to compare the error
performance of SuM-OFDM-IM with reference schemes. All simulations have been
performed under frequency-selective Rayleigh fading channels and perfect
channel estimation is assumed at the receiver side. OFDM parameters are
adjusted to $N=128$, $v=10$ and $L=16$ for all simulations. We define ”OFDM-
IM$(n,k)$ and CI-OFDM-IM$(n,k)$” as $k$ out of $n$ subcarriers being active in
a subblock, ”DM-OFDM$(n,k)$” as $k$ subcarriers out of $n$ exploiting the
primary $M$-ary PSK/QAM constellation, ”MM-OFDM-IM$(Q,M)$, SuM-OFDM-IM$(Q,M)$,
and S-SuM-OFDM-IM$(Q,M)$” with $M$ modes, each including $Q$ symbols.

As shown in Fig. 4, the theoretical BER curve obtained by (19) is an accurate
upper bound for simulation results of SuM-OFDM-IM for varying parameters
$(Q,M)$. Since (18) is an approximation, as spectral efficiency increases, the
number of terms in (19) also increases; therefore, the gap between simulation
and theoretical curves does not diminish.

Figure 5: Error performance comparison of SuM-OFDM-IM and S-SuM-OFDM-IM with
reference systems.

It is observed from Fig. 5 that SuM-OFDM-IM and S-SuM-OFDM-IM have superior
BER performances over reference systems for approximately the same spectral
efficiency. As seen from Fig. 5, SuM-OFDM-IM and S-SuM-OFDM-IM outperform
classical OFDM, OFDM-IM, DM-OFDM and MM-OFDM-IM. Additionally, at a BER value
of $10^{-5}$, SuM-OFDM-IM provides almost $10$ dB gain over MM-OFDM-IM
although providing a $12.5\%$ higher spectral efficiency. It is also shown
that LLR-based reduced ML detector provides the same error performance as that
of the ML detector. Although SuM-OFDM-IM performs only slightly worse in terms
of error performance than S-SuM-OFDM-IM, its spectral efficiency is $12.5\%$
higher.

Figure 6: Error performance comparison of SuM-OFDM-IM with CI-OFDM-IM.

In Fig. 6, SuM-OFDM-IM is compared with CI-OFDM-IM, which also provides a
second order diversity gain. First, we consider $2$ and $2.25$ bps/Hz for CI-
OFDM-IM and SuM-OFDM-IM, respectively. SuM-OFDM-IM provides a marginal gain at
a BER value of $10^{-5}$ than CI-OFDM-IM despite a $12.5\%$ improvement in
terms of spectal efficiency. Moreover, we consider $3$ and $3.25$ bps/Hz for
CI-OFDM-IM and SuM-OFDM-IM, respectively. SuM-OFDM-IM is able to provide a
noticeable gain at a BER value of $10^{-5}$ than CI-OFDM-IM despite yielding a
$8.33\%$ higher spectral efficiency. It can be clearly concluded that the
performance difference between SuM-OFDM-IM and CI-OFDM-IM increases with
spectral efficiency.

Figure 7: Error performance of SuM-OFDM-IM for different spectral efficiencies
with LLR detector for $n=4$.

In Fig. 7, the error performance of SuM-OFDM-IM with varying $(Q,M)$
parameters is given for different spectral efficiency values and $n=4$. LLR-
based reduced complexity ML detector is employed in all cases. As seen from
Fig. 7, the slopes of the curves are equal and the error performance degrades
with increasing spectral efficiency.

Figure 8: Error performance of SuM-OFDM-IM for different spectral efficiencies
with LLR detector for $n=8$.

Error performance of SuM-OFDM-IM for varying parameters $(Q,M)$ is shown in
Fig. 8 where $n=8$. We observe that the error performance is inversely
proportional to the size of the constellation. As $Q$ increases, the
percentage of higher rank values increases and it can be seen from Fig. 8 that
at a BER value of $10^{-5}$, the performances of schemes with $2$ bps/Hz and
$3$ bps/Hz are close to each other. By comparing Figs. 7 and 8 to investigate
the effect of $n$ on the BER performance, it can readily be deduced that the
system with $n=8$ and $3$ bps/Hz outperforms the system with $n=4$ and $3.25$
bps/Hz. This is because $r$ may even be equal to $8$ when $n=8$, resulting in
an improvement in error performance; nevertheless, the complexity increases
proportional to $n$ as seen in Table III. Additionally, it is more challenging
to reach a higher spectral efficiency for $n=8$ compared to $n=4$ by (7).
Therefore, we introduce interesting trade-offs among spectral efficiency,
decoding complexity and error performance.

Figure 9: Uncoded/Coded performance of OFDM, MM-OFDM-IM, SuM-OFDM-IM, and
S-SuM-OFDM-IM.

In Fig. 9, the error performance curves of uncoded and coded OFDM, MM-OFDM-IM,
SuM-OFDM-IM, and S-SuM-OFDM-IM are given. For all coded schemes, the
rate-$1/2$ convolutional code (CC) with octal generator sequence of [1, 3] and
a bit-interleaver are employed. As seen from Fig. 9, SuM-OFDM-IM and S-SuM-
OFDM-IM outperform OFDM at a BER value of $10^{-5}$ also with channel coding.
Moreover, it is observed that the difference between coded SuM-OFDM-IM and
coded OFDM increases for higher SNR values due to different slopes of their
curves. While the advantage of SuM-OFDM-IM is not significant compared to MM-
OFDM-IM for this specific case, SuM-OFDM-IM can also be a candidate for
potential future uncoded applications due to its diversity gain. However, in
the presence of coding, S-SuM-OFDM-IM outperforms MM-OFDM-IM. Therefore, we
present an interesting trade-off between error performance and spectral
efficiency for SuM-OFDM-IM and S-SuM-OFDM-IM.

TABLE III: Decoding Complexity Comparison System |  Detector |  Complexity Order
---|---|---
OFDM |  ML |  $\mathcal{O}(Q)$
OFDM-IM [4] |  Red. Comp. ML |  $\mathcal{O}(Q)$
DM-OFDM [19] |  Suboptimal |  $\mathcal{O}(Q_{A}+Q_{B})$
MM-OFDM-IM [23] |  Suboptimal |  $\mathcal{O}(\frac{Qn}{2}+\frac{Q}{2})$
CI-OFDM-IM [6] |  Red. Comp. ML |  $\mathcal{O}(Q^{2})$
SuM-OFDM-IM |  Red. Comp. ML |  $\mathcal{O}(\binom{n}{n/2}\frac{QM}{2})$
S-SuM-OFDM-IM |  Red. Comp. ML |  ${\mathcal{O}(\binom{n}{n/2}\frac{QM}{2})}$
TABLE IV: A Numerical Example For Complexity Comparison System |  Detector |  CMs per subcarrier
---|---|---
OFDM |  ML |  $4$
OFDM-IM |  Red. Comp. ML |  $4$
DM-OFDM |  Suboptimal |  $4$
MM-OFDM-IM |  Suboptimal |  $5$
CI-OFDM-IM |  Red. Comp. ML |  $64$
SuM-OFDM-IM |  Red. Comp. ML |  $48$
SuM-OFDM-IM |  ML |  $128$

###  V-A Complexity Analysis

In this subsection, we analyze the decoding complexity order of SuM-OFDM-IM
and make comparisons with the reference systems. A comparison for complexity
per subcarrier in terms of complex multiplications (CMs), is given in Table
III, where $Q_{A}$ and $Q_{B}$ stand for the primary and secondary
constellations for DM-OFDM, respectively, and $k$ is the number of active
subcarriers in a CI-OFDM-IM subblock. Due to space limitation, we do not give
the design of LLR-based reduced complexity ML detection for S-SuM-OFDM-IM,
however, it can be designed as in Algorithm 1 in a similar manner. In this
case, the complexity order of SuM-OFDM-IM and S-SuM-OFDM-IM would be
equivalent as seen Table III.

To further elaborate on the complexity order of the proposed system and
reference systems, we give a numerical example in Table IV. For a fair
comparison, the parameters in Figs. 6 and 5 are considered for CI-OFDM-IM with
$2$ bps/Hz and other systems, respectively. As seen from Table IV, we provide
an interesting trade-off among complexity and error performance. With these
specific parameters, although SuM-OFDM-IM is the second most complex system,
it provides a superior error performance as seen in Figs. 5 and 6. It should
also be noted that the LLR-based reduced complexity ML detector dramatically
reduces the complexity of SuM-OFDM-IM when compared to the ML detector.
Moreover, the LLR-based and low-complexity ML detector is not significantly
complex for varying $n$, $Q$ and $M$ values. For example, for a spectral
efficiency of $1.75$ bps/Hz, only $24$ CMs per subcarrier are required. Hence,
SuM-OFDM-IM is applicable in practice with the proposed low-complexity
detector.

Figure 10: Experimental testbed with USRP modules.

##  VI Experimental Results

In this section, we evaluate the performance of SuM-OFDM-IM in a practical
setup and give the models of the transmitter and receiver testbeds. As seen
from Fig. 10, we use two USRP 2943 modules for transmitter and receiver nodes
and an OctoClock CDA-2990 that generates $10$ MHz clock signal for
synchronization. For programming, virtual instrument (VI) components of
LabVIEW are utilized.

###  VI-A Transmitter Details

####  VI-A1 Channel Estimation

Comb-type channel estimation method with one dimensional linear interpolation
is employed to estimate the channel as in [22]. A certain number of
subcarriers are assigned as pilot tones. In this experiment, we consider the
ratio of data to pilot subcarriers as $4$, i.e, for $4$ data subcarriers, one
pilot tone is employed.

Figure 11: Overall symbol structure.

####  VI-A2 Symbol Structure

TABLE V: Experimental Parameters Carrier Frequency |  $2.45$ GHz
---|---
I/Q data rate |  $5\times 10^{5}$ samples/sec
Sampling rate $(f_{s})$ |  $6.25\times 10^{5}$ samples/sec
Number of data subcarriers |  $125$
Number of pilot tones |  $25$
Zero padding/FFT size $(N)$ |  $106/256$
CP size $(L)$ |  $32$
Distance between terminals |  $225$ cm

Overall symbol structure is shown in Fig. 11. We insert two identical training
sequences with $128$ samples for both TO and CFO estimation at the beginning
of the symbol. In the OFDM frame, the initial $53$, the center and the final
$52$ subcarriers are set to $0$. A total $150$ subcarriers are assigned to
data and pilot tones. Since the ratio of data subcarriers to pilot tones is
$4$, we have a total number of $125$ data and $25$ pilot subcarriers. After
the IFFT operation, a CP with $32$ samples is inserted in the time domain.
Experiment parameters are summarized in Table V. The bandwidth, subcarrier
spacing and symbol duration are $500$ kHz, $3.333$ kHz and $0.3$ ms,
respectively.

###  VI-B Receiver Details

In the presence of CFO $(\varepsilon)$ and TO $(\gamma)$, the received
baseband signal in time domain can be given as:

$\displaystyle{y}_{T}(\tau)=\frac{1}{N}\sum_{\beta=0}^{N-1}s_{F}(\beta)c_{F}(\beta)e^{j2\pi(\beta+\varepsilon)(\tau+\gamma)/N}+w_{T}(\tau),$
(23)

where $w_{T}(\tau)$, $\tau=1,2,\cdots,N$ represent the noise samples in the
time domain, whose vector form is given as
$\mathbf{w}_{T}=\mathrm{IFFT}\\{\mathbf{w}_{F}\\}$.

####  VI-B1 CFO and TO Estimation

After receiving signals, by using consecutive and identical training
sequences, timing offset (TO) and carrier frequency offset (CFO) are estimated
in time and frequency domain, respectively. By using sliding windows, which
have the same lengths as the training sequence, the similarity between two
blocks is obtained by autocorrelation [26]. Hence, TO can be ML estimated as:

$\hat{\gamma}=\underset{\gamma}{\arg\>\max}{\Bigg{\\{}\frac{\Big{|}\sum_{i=\gamma}^{N/2-1+\gamma}y_{T}(\tau+i)y_{T}^{*}(\tau+N/2+i)\Big{|}^{2}}{\Big{|}\sum_{i=\gamma}^{N/2-1+\gamma}y_{T}(\tau+N/2+i)\Big{|}^{2}}\Bigg{\\}}}.$
(24)

After obtaining the start of the frame (TO), by taking FFT of training
sequences, CFO is estimated in the frequency domain with the technique
proposed by Moose [27] as follows:

${\hat{\varepsilon}=\frac{1}{2\pi}\mathrm{tan}^{-1}\Bigg{\\{}\frac{\sum_{\beta=0}^{N-1}\Im\\{T_{1}^{*}(\beta)T_{2}(\beta)\\}}{\sum_{\beta=0}^{N-1}\Re\\{T_{1}^{*}(\beta)T_{2}(\beta)\\}}\Bigg{\\}}},$
(25)

where $T_{l},l=1,2$ is the $l$th received training sequence in the frequency
domain. Finally, the CP is removed and the CFO is compensated.

####  VI-B2 Channel Estimation

After performing FFT, data and pilot tones are separated. By employing pilot
tones, channel coefficients are estimated with one-dimensional linear
interpolation as in [22]. Finally, the estimated channel coefficients are fed
into the detector.

####  VI-B3 SNR Estimation

For SNR estimation, the method proposed in [28] is considered. The estimated
SNR values for all target systems can be given as

${\mathrm{SNR}}=\frac{\sum_{\beta=1}^{N}|y_{F}(\beta)|^{2}}{\sum_{\beta=1}^{N}|y_{F}(\beta)-s_{F}(\beta)h_{F}(\beta)|^{2}}.$
(26)

Figure 12: Error performance of OFDM, OFDM-IM and SuM-OFDM-IM under a
practical setup and with computer simulation.

###  VI-C Error Performance

In this subsection, we compare the error performance of SuM-OFDM-IM with OFDM
and OFDM-IM under a practical setup. ML detector is employed for all schemes.
We consider a three-tap Rician channel with uniform power delay profile and
$K$-factor of $10$ dB to compare simulation results with experimental results
similar to [22]. As seen from Fig. 12, SuM-OFDM-IM outperforms OFDM and OFDM-
IM in terms of error performance even though providing a $12.5\%$ higher
spectral efficiency. These experimental results are also consistent with the
simulation results.

##  VII Conclusion

In this paper, we have proposed SuM-OFDM-IM, which jointly determines MAPs and
SAPs by IM to convey additional information bits and replicates the
conventional QAM/PSK symbols to obtain a diversity gain. To enrich the
contribution, the separate selection of MAPs and SAPs have also been
investigated. A reduced complexity ML detector, which achieves the same error
performance as the ML detector, has been designed and a complexity analysis
has been given. Finally, the error performance of the proposed scheme has been
investigated both with computer simulations and SDRs under practical
impairments. Our exhaustive simulation and practical results demonstrate that
the proposed scheme remarkably outperforms OFDM and other OFDM-IM based
reference systems in terms of error performance while providing a high
spectral efficiency. Consequently, SuM-OFDM-IM can be a possible waveform
candidate for 5G and beyond eMBB services due to its appealing advantages such
as high spectral efficiency, reliability and diversity gain [29]. Other
precoding techniques can be employed to further increase the diversity order
of SuM-OFDM-IM. Moreover, the spectral efficiency can be further enhanced by
selecting more than two modes; however, it is left as a future study.

## References

  * [1] E. Basar, “Index modulation techniques for 5G wireless networks,” _IEEE Commun. Mag._ , vol. 54, no. 7, pp. 168–175, Jul. 2016.
  * [2] E. Basar, M. Wen, R. Mesleh, M. Di Renzo, Y. Xiao, and H. Haas, “Index modulation techniques for next-generation wireless networks,” _IEEE Access_ , vol. 5, pp. 16 693–16 746, 2017.
  * [3] T. Mao, Q. Wang, Z. Wang, and S. Chen, “Novel index modulation techniques: A survey,” _IEEE Commun. Surveys Tuts_ , vol. 21, no. 1, pp. 315–348, First Quarter 2019.
  * [4] E. Basar, Ü. Aygölü, E. Panayırcı, and H. V. Poor, “Orthogonal frequency division multiplexing with index modulation,” _IEEE Trans. Signal Process._ , vol. 61, no. 22, pp. 5536–5549, Nov. 2013\.
  * [5] Y. Xiao, S. Wang, L. Dan, X. Lei, P. Yang, and W. Xiang, “OFDM with interleaved subcarrier-index modulation,” _IEEE Commun. Lett._ , vol. 18, no. 8, pp. 1447–1450, Aug. 2014.
  * [6] E. Basar, “OFDM with index modulation using coordinate interleaving,” _IEEE Wireless Commun. Lett._ , vol. 4, no. 4, pp. 381–384, Aug. 2015.
  * [7] A. M. Jaradat, J. M. Hamamreh, and H. Arslan, “OFDM with subcarrier number modulation,” _IEEE Wireless Commun. Lett._ , vol. 7, no. 6, pp. 914–917, Dec. 2018.
  * [8] S. A. Nambi and K. Giridhar, “Lower order modulation aided BER reduction in OFDM with index modulation,” _IEEE Commun. Lett._ , vol. 22, no. 8, pp. 1596–1599, Aug. 2018.
  * [9] Y. Shi, K. Gao, J. Zhu, X. Lu, and S. Wang, “Orthogonal frequency division multiplexing with joint subblocks index modulation,” _IEEE Access_ , vol. 7, pp. 23 930–23 939, 2019.
  * [10] S. Dang, J. P. Coon, and G. Chen, “Adaptive OFDM with index modulation for two-hop relay-assisted networks,” _IEEE Trans. Wireless Commun._ , vol. 17, no. 3, pp. 1923–1936, Mar. 2018.
  * [11] S. Dang, G. Chen, and J. P. Coon, “Outage performance of two-hop OFDM with index modulation and multi-carrier relay selections,” _IEEE Wireless Commun. Lett._ , vol. 7, no. 6, pp. 926–929, Dec. 2018.
  * [12] E. Basar, “On multiple-input multiple-output OFDM with index modulation for next generation wireless networks,” _IEEE Trans. Signal Process._ , vol. 64, no. 15, pp. 3868–3878, Aug. 2016.
  * [13] B. Zheng, M. Wen, E. Basar, and F. Chen, “Multiple-input multiple-output OFDM with index modulation: Low-complexity detector design,” _IEEE Trans. Signal Process._ , vol. 65, no. 11, pp. 2758–2772, Jun. 2017.
  * [14] E. Öztürk, E. Basar, and H. A. Çırpan, “Multiple-input multiple-output generalized frequency division multiplexing with index modulation,” _Physical Communication_ , vol. 34, pp. 27–37, 2019.
  * [15] S. Lu, I. A. Hemadeh, M. El-Hajjar, and L. Hanzo, “Compressed-sensing-aided space-time frequency index modulation,” _IEEE Trans. Veh. Techol._ , vol. 67, no. 7, pp. 6259–6271, Jul. 2018.
  * [16] Q. Ma, P. Yang, L. Dan, X. He, Y. Xiao, and S. Li, “OFDM-IM-aided cooperative relaying protocol for cognitive radio networks,” in _2017 IEEE 18th Int. Workshop on Signal Processing Advances in Wireless Commun. (SPAWC)_ , July 2017, pp. 1–5.
  * [17] Q. Li, M. Wen, S. Dang, E. Basar, H. V. Poor, and F. Chen, “Opportunistic spectrum sharing based on OFDM with index modulation,” _IEEE Trans. Wireless Commun._ , pp. 1–1, 2019.
  * [18] E. Başar and E. Panayırcı, “Optical OFDM with index modulation for visible light communications,” in _2015 4th International Workshop on Optical Wireless Communications (IWOW)_ , Sep. 2015, pp. 11–15.
  * [19] T. Mao, Z. Wang, Q. Wang, S. Chen, and L. Hanzo, “Dual-mode index modulation aided OFDM,” _IEEE Access_ , vol. 5, pp. 50–60, Feb. 2017.
  * [20] J. Zheng and R. Chen, “Achieving transmit diversity in OFDM-IM by utilizing multiple signal constellations,” _IEEE Access_ , vol. 5, pp. 8978–8988, June 2017.
  * [21] T. Mao, Q. Wang, and Z. Wang, “Generalized dual-mode index modulation aided OFDM,” _IEEE Commun. Lett._ , vol. 21, no. 4, pp. 761–764, April 2017.
  * [22] S. Gokceli, E. Basar, M. Wen, and G. K. Kurt, “Practical implementation of index modulation-based waveforms,” _IEEE Access_ , vol. 5, pp. 25 463–25 473, 2017.
  * [23] M. Wen, E. Basar, Q. Li, B. Zheng, and M. Zhang, “Multiple-mode orthogonal frequency division multiplexing with index modulation,” _IEEE Trans. on Commun._ , vol. 65, no. 9, pp. 3892–3906, Sep. 2017.
  * [24] M. Wen, Q. Li, E. Basar, and W. Zhang, “Generalized multiple-mode OFDM with index modulation,” _IEEE Trans. Wireless Commun._ , vol. 17, no. 10, pp. 6531–6543, Oct. 2018.
  * [25] Q. Li, M. Wen, E. Basar, H. V. Poor, B. Zheng, and F. Chen, “Diversity enhancing multiple-mode OFDM with index modulation,” _IEEE Trans. Commun._ , vol. 66, no. 8, pp. 3653–3666, Aug. 2018.
  * [26] Y. S. Cho, J. Kim, W. Y. Yang, and C. G. Kang, _MIMO-OFDM Wireless communications with MATLAB_. John Wiley & Sons, 2010.
  * [27] P. H. Moose, “A technique for orthogonal frequency division multiplexing frequency offset correction,” _IEEE Trans. Commun._ , vol. 42, no. 10, pp. 2908–2914, Oct 1994.
  * [28] H. Arslan and S. Reddy, “Noise power and SNR estimation for OFDM based wireless communication systems,” in _Proc. of 3rd IASTED Int. Conf. on Wireless Optical Commun. (WOC), Banff, Alberta, Canada_ , 2003.
  * [29] S. Dang, O. Amin, B. Shihada, and M.-S. Alouini, “What should 6G be?” _Nature Electronics_ , vol. 3, no. 1, pp. 20–29, 2020.

|  Ali Tugberk Dogukan (S’20) received the B.S. degree from Istanbul Technical
University, Istanbul, Turkey, in 2018, and the M.S. degree from Koç
University, Istanbul, Turkey, in 2020. He is currently a Research and Teaching
Assistant at Koç University while pursuing his Ph.D. degree at the same
university. His research interests include waveform design, signal processing
for communications, index modulation, and software defined radio-based
practical implementation. He was a Reviewer for IEEE.
---|---
|  Ertugrul Basar (S’09-M’13-SM’16) received the B.S. degree (Hons.) from
Istanbul University, Turkey, in 2007, and the M.S. and Ph.D. degrees from
Istanbul Technical University, Turkey, in 2009 and 2013, respectively. He is
currently an Associate Professor with the Department of Electrical and
Electronics Engineering, Koç University, Istanbul, Turkey and the director of
Communications Research and Innovation Laboratory (CoreLab). His primary
research interests include MIMO systems, index modulation, intelligent
surfaces, waveform design, visible light communications, and signal processing
for communications. Recent recognition of his research includes the IEEE
Communications Society Best Young Researcher Award for the Europe, Middle
East, and Africa Region in 2020, Science Academy (Turkey) Young Scientists
(BAGEP) Award in 2018, Turkish Academy of Sciences Outstanding Young Scientist
(TUBA-GEBIP) Award in 2017, and the first-ever IEEE Turkey Research
Encouragement Award in 2017. Dr. Basar currently serves as a Senior Editor of
the IEEE Communications Letters and the Editor of the IEEE Transactions on
Communications, Physical Communication (Elsevier), and Frontiers in
Communications and Networks.
---|---


# Accelerating Airy tensor modes of cosmological gravitational waves

Claudio Aravena-Plaza<EMAIL_ADDRESS>(corresponding author)
Departamento de Física, Facultad de Ciencias, Universidad de Chile, Santiago
7800003, Chile.  Víctor Muñoz  Departamento de Física, Facultad de Ciencias,
Universidad de Chile, Santiago 7800003, Chile.  Felipe A. Asenjo  Facultad de
Ingeniería y Ciencias, Universidad Adolfo Ibáñez, Santiago 7491169, Chile.

###### Abstract

From a classical analysis, it is shown that the nondiffractive accelerating
gravitational Airy wave packets are solutions of Einstein equations for their
linearized tensor modes in a Friedmann-Lemaître-Robertson-Walker cosmological
background filled with a perfect fluid, with equations of state $w=1/3$ and
$w=-1/3$. These solutions have finite energy, presenting accelerating behavior
due to the structured spatial form of the wavepacket. This is manifested by
curved trajectories along the wave path. Also, using spectral functions, it is
possible, with these packets, to construct more general, arbitrary wave
packets. All these new solutions bring insights on new forms for gravitational
wave propagation.

##  I Introduction

Nondiffractive wave packets are known solutions of the Schrödinger equation
berry ; Yuce ; Unnikrishnan  and Maxwell’s equations kaminer ; Hacyan .
Several experimental studies have been performed to show them as physical
modes, for instance, for light propagation kaminer ; Hacyan ; Aleahmad ;
Chremmos ; Efremidis ; Patsyk ; Bandres ; Esat , fluids shenhe , sound
zhao2014delivering ; chen2019broadband ; tang2021generation , heat diffusion
functions2004applications ; asenjo2021accelerating ; asenjo2022airy  and
relativistic electron plasmas winkler2023exact . In gravitational waves, these
kinds of solutions have also been obtained in flat spacetimes Asenjo . The
solutions in terms of Airy modes represent accelerated wavepackets without
diffraction in their trajectory, due to a non-vanishing Bohm potential
hojman2020classical , and thus, they describe solutions different to those
described through plane waves.

On the other hand, from the linearization of the Einstein equations for a
Friedmann-Lemaître-Robertson-Walker (FLRW) cosmological background metric, the
dynamical equations for the scalar, vector and tensor modes of gravitational
perturbations are straightforwardly obtained misner . The tensor modes are of
particular importance because they describe gravitational waves on this
cosmology. In the simplest case, those tensor modes propagate in a Universe
filled with a perfect fluid with a specific equation of state. Depending on
that, solutions for different stages of the Universe are modeled weinberg .
The detection of these waves is important due to the idea of an inflationary
early Universe producing primordial gravitational waves osti_981884 .

In general, the study of tensor modes for an FLRW background metric can be
performed exactly. Their dynamical evolution is determined by the equation of
state of the cosmological content of the Universe, mainly described by a
perfect fluid with pressure $p$ and its energy density $\rho$, and thus, an
equation of state $w=p/\rho$. The different values of $w$ characterize the
evolution of the scale factor of the Universe. The distinction is usually made
between oscillatory solutions with a cutoff value of $w>-1/3$, where the
solutions oscillate and damp, and with $w<-1/3$, where the oscillation occurs
and then freezes at some value.

In this work, we focus in cosmologies with specific values $w=1/3$ and
$w=-1/3$. In the following, we show that in these two particular cosmological
scenarios, the propagation of gravitational waves can be described in terms of
Airy functions. This behavior is due only to the structured spatial form of
the wave, transverse to its propagation along cosmological light-cone
coordinates. This produces accelerating curved (parabolic) trajectories for
the gravitational wave propagation. Both values of $w$ here considered define
different forms in which the Universe evolves. The characterization of a
Universe with $w=1/3$ represents a stage of the Universe where
ultrarelativistic particles predominate. In this case, primordial
gravitational waves can be found when Einstein linearized equations are
coupled with Vlasov’s equations khlebnikov1984long . On the other hand, a
Universe with $w=-1/3$ falls within the so-called dark energy cosmologies, and
consists of linear-type forms of energy, called the cosmic string. A cosmic
string could arise as a topological defect that occurred during an early phase
transition in the Universe. The NanoGrav Collaboration has recently reported
strong evidence of common spectrum process, which was interpreted as a
stochastic gravitational wave background (SGWB) in the framework of cosmic
string PhysRevD.107.042003 . The importance of this measurement opens the
field of observation of primordial background waves nanograv . Furthermore,
cosmic strings can produce primordial gravitational waves, such as
gravitational wave bursts and SGWB Siemens_2007 ;
https://doi.org/10.48550/arxiv.2207.03510 ;
https://doi.org/10.48550/arxiv.2006.13872 ; CMB ; Emond_2022 ; Caprini_2018 .
In general, to show that a Universe with a cosmic string can produce
primordial gravitational waves, different models described by some action are
used, for example, using the Nambu-Goto action auclair2021window  or
Lagrangian models Emond_2022 . Unlike those works, here we perform a tensor
modal analysis finding Airy-type wave solutions.

The studied gravitational wave modes in this work correspond to spacetime
perturbations of the FLRW background metric, which is described by the
spacetime interval

$\displaystyle
ds^{2}=g_{\mu\nu}dx^{\mu}dx^{\nu}=a^{2}\left(-d\eta^{2}+\gamma_{ij}dx^{i}dx^{j}\right)\,,$
(1)

where $a=a(t)$ is the scale factor of the Universe, and $x^{\mu}=(t,x^{i})$
(with $\mu,\nu=0,1,2,3$, and $i,j=1,2,3$). We have chosen the speed of light
$c=1$. We can describe the dynamics in terms of the conformal cosmological
time $\eta=\int dt/a$. The three-dimensional spatial metric $\gamma_{ij}$ can
be represented in spherical coordinates as
$\gamma_{ij}dx^{i}dx^{j}=dD^{2}+D_{A}^{2}d\Omega$, where $D$ is the comoving
distance, $D_{A}=K^{-1/2}\sin(K^{1/2}D)$ is the angular diameter distance, and
$d\Omega$ is the differential angular element. $K$ is a constant representing
the spatial curvature of the Universe ($K=0$ for spatially flat, $K=1$ for
closed Universe, and $K=-1$ for open Universe).

From the above metric, the gravitational tensor modes $H_{T}$ can be
straightforwardly obtained weinberg . They are obtained as a perturbation of
the spatial components of metric (1). By defining the effective gravitational
wave tensor mode $Y$ (for any component) by the relation $H_{T}=Y/a$, we get
the following equation for the cosmological gravitational tensor modes (see
Appendix):

$-\nabla^{2}Y+\ddot{Y}+\left(2K-\frac{\ddot{a}}{a}\right)Y=0\,,$ (2)

where the overdot stands for the conformal time derivative, $d/d\eta$. Also,
the tensor $Y$ must satisfy a gauge condition
$\partial_{\beta}Y^{\beta}_{~{}\alpha}=0$.

Furthermore, the background FLRW metric evolves following the equation
${\ddot{a}}/{a}=({4\pi G}/{3})a^{2}\rho(1-3w)-K$, where $G$ is the
gravitational constant. The fluid content satisfies the continuity equation
$\dot{\rho}/\rho=-3\left(1+w\right)\left(\dot{a}/a\right)$, whose solution is
$\rho=\rho_{0}a^{-3(1+w)}$, where $\rho_{0}$ is the initial energy density.
Substituting all the above in Eq. (2) for a spatially flat Universe with
$K=0$, we finally obtain

$-\nabla^{2}Y+\ddot{Y}-\frac{4\pi G}{3}a^{-1-3w}\rho_{0}(1-3w)Y=0\,.$ (3)

The solution of this equation can be found exactly in terms of Bessel
functions for any equation of state different to $-1/3$ and $1/3$.

##  II Airy cosmological gravitational waves

The main purpose of this work is to show that, from Eq. (3), Airy-type
wavepacket solutions are obtained when $w=1/3$ and $w=-1/3$. This is performed
in the following sections, showing and discussing the different possible
solutions. In both cases, the wave propagates along curved trajectories in the
plane formed by a cosmological conformal light-cone coordinate and a
transverse direction.

###  II.1 Accelerating gravitational waves for $w=1/3$

For the case of the equation of state for ultra relativistic matter, wave
equation (3) reduces to

$-\nabla^{2}Y+\ddot{Y}=0\ .$ (4)

It is customary to solve this equation in terms of a sinusoidal basis. By
considering a polarized wave given by $Y=Y_{zz}=Y_{zz}(\eta,x,y)$, which
satisfies the Lorenz gauge condition $\partial_{z}Y_{zz}=0$, we can solve (4)
as a plane wave $Y=\exp(i\omega\eta-ik_{x}x-ik_{y}y)$, such that
$\omega=\sqrt{k_{x}^{2}+k_{y}^{2}}$, implying a light-cone propagation of the
wave. This is achieved by the symmetric aspects of spatial coordinates in this
solution. However, other solutions can be constructed where the different
spatial coordinates are not treated in the same footing. These solutions lead
to an accelerated propagation of the wave, whose amplitude can be written in
terms of Airy functions. For this, we assume that the polarized wave has the
form

$Y=Y_{zz}(\eta,x,y)=\Upsilon(\zeta,y)\exp(ik\tau),$ (5)

where $k$ is an arbitrary constant playing the role of a wavenumber. Also, we
have defined the light-cone coordinates

$\displaystyle\zeta=x-\eta\,,\qquad\tau=x+\eta\,.$ (6)

Using this ansatz in Eq. (4), the equation for $\Upsilon$ is

$4ik\frac{\partial\Upsilon}{\partial\zeta}+\frac{\partial^{2}\Upsilon}{\partial^{2}y}=0\,,$
(7)

which is equivalent to a Schrödinger equation for a free particle. It is very
well-known that a possible solution of this equation can be written in terms
of Airy functions Ai berry , being

$\Upsilon(\zeta,y)=\text{Ai}\left(2ky-k^{2}\zeta^{2}\right)\exp\left(2ik^{2}y-\frac{2i}{3}k^{3}\zeta^{3}\right)\,.$
(8)

This solution represents a nondiffractive accelerating gravitational wave
packet, such that its intensity bends its trajectory in a parabolic curve in
the $\zeta$-$y$ plane.

However, solution (8) has infinite energy. A finite energy Airy gravitational
wavepacket has already been studied in Ref. (Asenjo, ). Thus, for this case,
the cosmological gravitational wave, with finite energy, is given by

$\Upsilon(\zeta,y)=\text{Ai}\left(2ky-k^{2}\zeta^{2}+ibk\zeta\right)\\\
\times\exp\left(bky-
bk^{2}\zeta^{2}+\frac{i}{4}b^{2}k\zeta+2ik^{2}y\zeta-\frac{2i}{3}k^{3}\zeta^{3}\right),$
(9)

in terms of an arbitrary real parameter $b$ that allows the normalization of
the solution. This wavepacket has finite energy, since $\int
dy|\Upsilon|^{2}=\int dy~{}\Upsilon^{*}\Upsilon$ is independent of $\zeta$.
Then, it is obtained that the wave packets are square integrable along
$\zeta=0$ $(x=\eta)$ as

$\int_{-\infty}^{\infty}dy|\Upsilon(\zeta,y)|^{2}=\int_{-\infty}^{\infty}dy|\Upsilon(0,y)|^{2}=\frac{e^{b^{3}/12}}{4k\sqrt{b\pi}},$
(10)

implying its finite energy.

In Fig. 1(a), we plot the magnitude of the solution (9), for $b=1$ and $k=1$.
The parabolic trajectory followed by the intensity of the wave appears in the
$\zeta$-$y$ plane. In red line, the maximum maximorum of the solution (9) is
highlighted. For this case, and for small $\zeta$ values, the parabolic
trajectory of this maximum is described by the approximate equation
$y\approx{\zeta^{2}}/{3}+{\zeta}/{40}-{1}/{5}$, such that the initial
acceleration of the wave, along this parabolic trajectory on the light-cone
coordinates, is given by

$\frac{d^{2}y}{d\zeta^{2}}\approx\frac{2}{3}\,.$ (11)

The behavior of the different parts of the gravitational wavepacket (9) is
non-local, making this solution very different to a plane wave. Instead of
following straight lines in an isotropic space, the different parts of this
solution (in particular, the maximum intensity of the gravitational wave) bend
its trajectory.

Figure 1: (a) Density plot for the magnitude of $\Upsilon$, given by Eq. (9),
for $b=1$ and $k=1$. In red line, the parabolic trajectory of the maximum
intensity of the solution is shown. (b) Density plot for the magnitude of
generalized wavepacket ${\cal Y}(\zeta,y)$, given in Eq. (17), for spectral
function $\varrho=\exp\left(-(k-k_{0})^{2}\right)$, with $b=1$ and $k_{0}=1$.

###  II.2 Accelerating gravitational waves for $w=-1/3$

For this case of the equation of state, the Eq. (3) becomes

$-\nabla^{2}Y+\ddot{Y}-\frac{8\pi G}{3}\rho_{0}Y=0.$ (12)

Thus, the content of the Universe acts as en effective imaginary mass for
gravitational propagation.

The Airy form for this propagation can be obtained using an ansatz similar to
the previous section. We consider a polarized wave $Y_{zz}$ with the form

$Y=Y_{zz}(\eta,x,y)=\Upsilon(\zeta,y)\exp\left[ik\tau+iS(\zeta)\right]\,,$
(13)

in terms of light-cone coordinates, and transverse coordinate $y$ (again, $k$
is an arbitrary constant). With this, Eq. (12) reduces to

$4ik\frac{\partial\Upsilon}{\partial\zeta}+\frac{\partial^{2}\Upsilon}{\partial^{2}y}-4\left(k\frac{dS}{d\zeta}-\frac{2\pi
G}{3}\rho_{0}\right)\Upsilon=0.$ (14)

In order to have Airy-type wave packet solutions, the function $S$ must be
chosen to be

$S(\zeta)=\frac{2\pi G\rho_{0}}{3k}\zeta\,,$ (15)

to finally get

$4ik\frac{\partial\Upsilon}{\partial\zeta}+\frac{\partial^{2}\Upsilon}{\partial^{2}y}=0.$
(16)

This equation is the same as the one obtained for the $w=1/3$ case. Therefore,
it has the same mathematical finite energy solution (9), evolving, for
instance, as it is depicted in Fig. (1)(a).

Therefore, the differences between the $w=-1/3$ and $w=1/3$ cases appear only
on the phases of the gravitational waves, and not in their accelerating curved
trajectories.

##  III Discussion on General wavepackets

The Airy non-diffractive wavepacket solutions has been studied and measured in
quantum mechanics, electromagnetism and optics berry ; Yuce ; kaminer ; Hacyan
; Aleahmad ; Efremidis ; Patsyk ; Bandres ; Esat ; Unnikrishnan , and the
solutions presented in this work are the extension of this natural phenomena
to cosmological gravitational waves. They are particularly possible in the two
distinctive Universes with $w=1/3$ and $w=-1/3$.

Importantly, the above solutions have a degree of freedom that can be used to
construct a more general set of wavepackets solutions of this kind. This can
be done by constructing the wavepacket as an average of the previous solution
through a spectral function that depends only on the arbitrary parameter $k$.
The general wavepacket acquires the form

${\cal Y}(\zeta,y)=\int dk\,\varrho(k)\,Y(\zeta,y,k)\,,$ (17)

where $Y(\zeta,y,k)$ is any of the solutions in the previous sections, and
$\varrho(k)$ is an arbitrary spectral function, that can be chosen at will to
define different wavepacket structures. Note that this wavepacket definition
is clearly a solution of the equations (4) and (12), for the corresponding
solutions for $Y$. Thus, this wavepacket correspond to an average on arbitrary
distributions of different wavenumbers $k$.

As an appropriate example to demonstrate the dynamics of this general
wavepacket (17), we plot its magnitude in Fig. (1)(b), for a spectral function
determined by the Gaussian form $\varrho(k)=\exp\left(-(k-k_{0})^{2}\right)$.
This particular spectral function selects only wavenumbers close to $k_{0}$ as
contributors to the Airy dynamics. As a consequence, the Airy wavepacket
becomes more Gaussian as $\zeta$ progresses, such that the accelerating
features of the wave persist almost only for the maximum intensity lobe. This
example shows how an almost-Gaussian accelerating gravitational wavepacket can
be constructed from a set of Airy wavepackets. In general, any desired form of
accelerating gravitational wavepacket can be obtained by using an appropriate
spectral function.

The two cosmological scenarios, for fluids with equations of state $w=1/3$ and
$-1/3$, are not commonly considered into the study of cosmological
gravitational waves. With this work, we hope to bring new insights on the
accelerating features of the gravitational wavepackets on those cosmological
settings.

###### Acknowledgements.

CAP thanks to ANID No. 7322/2020. This project has been financially supported
by FONDECyT under contracts No. 1201967 (VM), and No. 1230094 (FAA).

##  IV APPENDIX

The gravitational tensor modes $H_{T}$ for a spacetime perturbation of the
FLRW metric follow the wave equation weinberg

$\nabla^{2}H_{T}-{a^{2}}\left(\dfrac{d^{2}}{dt^{2}}H_{T}\right)-3{a}\left(\dfrac{d}{dt}a\right)\left(\dfrac{d}{dt}H_{T}\right)=0\,.$
(18)

The evolution of these modes as a function of the cosmological time $\eta$ can
be obtained by noticing that $\dot{a}\equiv{da}/{d\eta}=a({da}/{dt})$.
Incorporating curvature constant $K$, we obtain (where
$\dot{H}_{T}\equiv{d}H_{T}/d\eta$)

$-\nabla^{2}H_{T}+\ddot{H}_{T}+2\frac{\dot{a}}{a}\dot{H}_{T}+2KH_{T}=0\,.$
(19)

## References

  * [1] N.L Balazs M.B. Berry.  Nonspreading wave packets.  Am. J. Phys., 47:264, 1979.
  * [2] C. Yuce.  Self-accelerating parabolic cylinder waves in 1-D.  Phys. Lett., 380, 2016.
  * [3] K. Unnikrishnan and A. R. P. Rau.  Uniqueness of the Airy packet in quantum mechanics.  American Journal of Physics, 64, 1996.
  * [4] I. Kaminer, R. Bekenstein, J. Neminovsky, and M. Segev.  Nondiffracting accelerating wave packets of Maxwell’s equations.  Phys. Rev. Lett., 108:163901, 2012.
  * [5] S. Hacyan.  Relativistic accelerating electromagnetic waves.  J. Opt, 13, 2011.
  * [6] P. Aleahmad, H. Moya Cessa, I. Kaminer, and M. Segev. D.N. Christodoulides.  Dynamics of accelerating Bessel solutions of maxwell’s equations.  J. Opt. Soc. Am. A, 33, 2016.
  * [7] N.K Efremidis I.D. Chremmos.  Nonparaxial accelerating Bessel-like beams.  Phys. Rev. A, 88:063816, 2013.
  * [8] N.K Efremidis, Z Chen, M. Segev, and D.N Christodoulies.  Airy beams and accelerating waves: an overview of recent advances.  Optica, 6, 2019.
  * [9] A. Patsyk, M.A. Bandres, R. Bekenstein, and M. Segev.  Observation of accelerating wave packets in curved space.  Phys. Rev X, 8:011001, 2018.
  * [10] J.C. Gutiérrez-Vega M.A Bandres.  Airy-Gauss beams and their transformation by paraxial optical systems.  Opt. Exp, 15, 2007.
  * [11] A.F. Abouraddy H. Esat Kondakci.  Airy wave packets accelerating in space-time.  Phys. Rev. Lett, 120, 2018.
  * [12] Shenhe Fu, Yuval Tsur, Jianying Zhou, Lev Shemer, and Ady Arie.  Propagation dynamics of Airy water-wave pulses.  Phys. Rev. Lett, 115, 2015.
  * [13] Sipei Zhao, Yuxiang Hu, Jing Lu, Xiaojun Qiu, Jianchun Cheng, and Ian Burnett.  Delivering sound energy along an arbitrary convex trajectory.  Scientific Reports, 4(1):6628, 2014.
  * [14] Di-Chao Chen, Xing-Feng Zhu, Da-Jian Wu, and Xiao-Jun Liu.  Broadband Airy-like beams by coded acoustic metasurfaces.  Applied Physics Letters, 114(5), 2019.
  * [15] Shuai Tang, Bin Ren, Yuxin Feng, Jie Song, and Yongyuan Jiang.  The generation of acoustic Airy beam with selective band based on binary metasurfaces: Customized on demand.  Applied Physics Letters, 119(7), 2021.
  * [16] Olivier Vallée and Manuel Soares.  Airy funcions. applications to physics, 2004.
  * [17] Felipe A Asenjo and Sergio A Hojman.  Accelerating solutions to diffusion equation.  The European Physical Journal Plus, 136:1–8, 2021.
  * [18] Felipe A Asenjo and Sergio A Hojman.  Airy heat bullets.  The European Physical Journal Plus, 137(10):1–9, 2022.
  * [19] Maricarmen A Winkler, Camilo Vásquez-Wilson, and Felipe A Asenjo.  Exact and paraxial airy propagation of relativistic electron plasma wavepackets.  The European Physical Journal D, 77(6):97, 2023.
  * [20] S.A Hojman F.A Asenjo.  Nondiffracting gravitational waves.  Eur. Phys. J.C, 81, 2021.
  * [21] Sergio A Hojman and Felipe A Asenjo.  Classical and quantum dispersion relations.  Physica Scripta, 95(8):085001, 2020.
  * [22] J.A. Wheeler C. W. Misner, K. S. Thorne.  Gravitation.  Princeton, 1973.
  * [23] Steven Weinberg.  Cosmology.  Oxford, 2008.
  * [24] Lawrence M. Krauss, Scott Dodelson, and Stephan Meyer.  Primordial gravitational waves and cosmology.  Science, 328, 2010.
  * [25] VI Khlebnikov.  Long gravitational waves in an ultrarelativistic two-component universe.  Zhurnal Eksperimental’noi i Teoreticheskoi Fiziki, 86, 1984.
  * [26] Y. M. Wu, Z. C. Chen, and Q. G. Huang.  Search for stochastic gravitational-wave background from massive gravity in the nanograv 12.5-year dataset.  Phys. Rev. D, 107, (2023).
  * [27] M. Benetti, G. Leila, and S. Vagnozzi.  Primordial gravitational waves from nanograv: A broken power-law approach.  Physical review D: Particles and fields, 105, (2021).
  * [28] X. Siemens, V. Mandic, and J. Creighton.  Gravitational-wave stochastic background from cosmic strings.  Physical review letters, 98(4), (2007).
  * [29] P. Auclair, S. Blasi, V. Brdar, and K. Schmitz.  Gravitational waves from current-carrying cosmic strings, (2022).
  * [30] R. Zhou and L. Bian.  Gravitational waves from cosmic strings and first-order phase transition, (2020).
  * [31] A. Bodas and R. Sundrum.  Large primordial fluctuations in gravitational waves from phase transitions, (2022).
  * [32] W. T. Emond, S. Ramazanov, and R. Samanta.  Gravitational waves from melting cosmic strings.  Journal of Cosmology and Astroparticle Physics, (2022), (2022).
  * [33] C. Caprini and D. G. Figueroa.  Cosmological backgrounds of gravitational waves.  Classical and Quantum Gravity, 35, (2018).
  * [34] P. Auclair, K. Leyde, and D. A. Steer.  A window for cosmic strings.  Journal of Cosmology and Astroparticle Physics, (2023), 2023.


# A potential system with infinitely many critical periods

Jihua Wang<EMAIL_ADDRESS>(J. Wang)

###### Abstract.

In this paper, we propose an analytical non-polynomial potential system which
has infinitely many critical periodic orbits in phase plane. By showing the
existence of infinitely many $2\pi-$ periodic solutions, the proof bases on
variational methods and the properties of Bessel function. The result provides
an affirmative example to Dumortier’s conjecture [Nonlinear Anal. 20 (1993)].

###### Key words and phrases:

Period function; Potential system; Critical point; Calculus of Variation.

###### 1991 Mathematics Subject Classification:

34C23, 34C25, 37G15

##  1\. Introduction and main result

The phenomenon of Oscillations (or vibrations) is popular in dynamical system.
The study on its period (or the frequency) is very important in the
applications of engineer or the other technology. Consider the conservative
Newtonian model $\ddot{x}+g(x)=0$ or equivalently in the phase plane

(1.1) $\begin{split}\frac{dx}{dt}=-y,\quad\frac{dy}{dt}=g(x).\end{split}$

It generates the classical one degree of freedom Hamiltonian vector fields on
the plane with a “kinetic + potential” Hamiltonian of the form
$H(x,y)=\frac{y^{2}}{2}+G(x)=h$, where $G(x)=\int_{0}^{x}g(s)ds$ stands for
the potential energy.

Assume that the origin $O$ is the center of system (1.1), the biggest
punctured neighborhood enclosing which is called the period annulus. In what
follows we denote it by $\mathcal{P}$ and the projection of which onto
$x-$axis by $(x_{l},x_{r})$. Note that the annulus $\mathcal{P}$ is foliated
by periodic orbits. Assigning to each periodic oval in $\mathcal{P}$ its
minimal period of the motion along it, then we define period function $T$ of
the center. Since the period function is defined on the set of periodic orbits
in $\mathcal{P}$, it is natural firstly to parameterize the set. Usually we
parameterize it by the energy level $h$ of the periodic oval and denote the
oval by $\Gamma_{h}$ accordingly. It is known that for a given continuum of
periodic orbits, the number, character (maximum or minimum) and distribution
of critical periods do not depend on the particular parametrization of the set
of periodic orbits applied. If the period function $T(h)$ is constant in a
neighborhood of the center, the center is called isochronous. The isolated
zeros $\hat{h}$ of the first order derivative of $T$ are called critical
periods, accordingly the periodic oval $\Gamma_{\hat{h}}$ is called critical
periodic orbit. The study on period function is closely related to two-point
nonlinear boundary-value problem of ODEs [4] and some bifurcation problems
such as subharmonic bifurcations of periodic oscillations and the bifurcation
of steady-state solutions of a reaction-diffusion equation [5, 6, 7].

As commented in [11, 12], many problems related to the period function and its
critical point are the counterpart of that about the poincaré return map and
limit cycle (isolated periodic orbit), which is the main focus of the
celebrated Hilbert’s 16th problem. In 1980, Z. Zhang [17] found an example
that a planar analytic system has an infinite number of limit cycles which
answered a conjecture opened for many years. Analytical property of vector
field only assure the local finiteness of number of limit cycles for a planar
system of differential equation while it could not deduce the global
finiteness of limit cycles.

Chicone and Dumortier [8] proved that the critical periodic orbits of an
analytic center can not accumulate on a compact set. If potential $g$ of
System (1.1) is nonlinear polynomial, then the number of critical periods is
finite globally. To the best of our knowledge, there are few papers on the
existence of multiple (no less than two) critical periods [10, 12, 13, 14]
which should be taken into account when looking for multiple solutions of
boundary value problems [9]. F. Dumortier remarked in [8] “it is possible that
an analytic vector field has an infinite number of critical periods. For
example, a linear vector field on the plane with a period annulus has a
constant period function”. While just we define the limit cycle which refers
to the isolated periodic orbit, in present paper we intend to propose an
nontrivial example which has an infinite number of critical periods
demonstrated as the strict extremum of period function.

Consider system (1.1) with potential $g(x)=x+\lambda\sin(x)$.

(1.2)
$\begin{split}\frac{dx}{dt}=-y,\quad\frac{dy}{dt}=x+\lambda\sin(x),\end{split}$

where the parameter $\lambda\in\,[-1,+\infty)$ is real and the potential
energy

$G(x)=\int_{0}^{x}\,g(s)ds=\frac{x^{2}}{2}+\lambda(1-\cos(x))$

is even. For all $\lambda\geq\,-1$ the period annulus is the whole plane,
accordingly the origin is called a global center. The periodic orbit is the
level curve
$\\{(x,y)\in\,R^{2}|\,\frac{y^{2}}{2}+G(x)=h,\;h\in\,(0,+\infty)\\}.$ When
$\lambda\in(-1,+\infty),$ the center is elementary at which the Jacobian
Matrix has a pair of pure imaginary eigenvalues, i.e., the linearized matrix
is non-degenerate. When $\lambda=-1$, the origin is a nilpotent center.

The system (1.1) with potential $g(x)=x+\lambda\sin(x)$ could be regarded as a
perturbation of the rigidly isochronous center at the origin (i.e.,
$\dot{\theta}\equiv\,1$ in the polar coordinates, as defined in [3]). In
another hand C. Chicone [4] consider the Neumann boundary value problem

(1.3) $\begin{split}\ddot{x}+\mu\,\sin x=0,\\\
\dot{x}(0)=0,\quad\dot{x}(T)=0.\end{split}$

with $\mu>0$. Here $x\equiv\,0$ is always a solution. One may be interested in
subharmonic solutions with $N$ nodes, that is $x(t)=0$ exactly $N$ times on
$(0,T)$. In geometrical viewpoint the period annulus enclosing the origin has
the minimum period $\frac{2\pi}{N}$. It reveals that the analysis of boundary
value problem can be completely understood in terms of period function which
assigns to each $x\in\,(0,\pi)$. It is closely related to the problem of
critical periods.

In present paper, we are interested in the bifurcation of critical periods,
roughly speaking, the emergence or disappearance of critical period as we
perturb a period annulus enclosing a center. One could get critical period
bifurcation by perturbing the periodic orbits of isochronous systems [12, 14],
we ask for the number of critical period orbits that persists. In a similar
manner as in studying limit cycle bifurcation whether if the period orbit
persists after a small perturbation of the continua of periodic orbits
enclosing a center which is known as the weakened Hilbert’s 16th problem first
posed by Arnold [2].

The main result is stated as follows

######  Theorem 1.1.

The origin of potential system (1.2) has infinitely many critical periods, it
appears alternatively as maximum and minimum.

##  2\. Preliminary lemmas and results

The asymptotic expansion of period function near the simple center is well
known, See Theorem 4.3 of [8] or Proposition 3.5 of [16] for instance.

######  Lemma 2.1.

If $\lambda\in(-1,+\infty)$, the origin $O(0,0)$ of potential system (1.2) is
a simple center, period function $T(h)$ can be analytically extended to the
origin, the period function and its first order derivation $T^{\prime}(h)$ are
stated as follows

(2.1)
$T(0^{+})=\frac{2\pi}{\sqrt{g^{\prime}(0)}},\;T^{\prime}(0^{+})=\pi\frac{{5(g^{\prime\prime}(x))^{2}-3\,g^{\prime}(x)g^{\prime\prime\prime}(x)}}{12(g^{\prime}(x))^{\frac{7}{2}}}|_{x=0}.$

If $\lambda=-1,$ the origin $O(0,0)$ is the nilpotent center of order 1, the
period function $T$ has a pole of order 1 at $O$, that is $T(0^{+})=+\infty$
and $T^{\prime}(0^{+})=-\infty.$

When the period annulus enclosing a center is the whole plane, the center is
global. Concerning the asymptotical behaviors of $T(h)$ near the outer
boundary of $\mathcal{P}$, we have

######  Lemma 2.2.

([15]) If potential $g$ of system (1.1) is locally Lipschitz and satisfies

  1. (1).

strict condition of signs: $xg(x)>0$ for all $x\neq\,0$;

  2. (2).

semilinear condition:
$0<\liminf\limits_{|x|\rightarrow\,+\infty}\frac{g(x)}{x}\leq\,\limsup\limits_{|x|\rightarrow\,+\infty}\frac{g(x)}{x}<+\infty,$

then the period function $T(h)$ satisfies

$0<\liminf\limits_{h\rightarrow\,+\infty}T(h)\leq\,\limsup\limits_{h\rightarrow\,+\infty}T(h)<+\infty.$

On studying the isochronicity of potential center, Urabe’s striking conclusion
[6] states as follows

######  Lemma 2.3.

If the potential $g$ is odd function, the only type of planar potential system
(1.1) having an isochronous center at the origin is the linear one.

Concerning necessary condition to the monotonicity of period for even
potential system (1.1), Opial assert in Theorem 6 in [5] it is stated as
follows.

######  Lemma 2.4.

If the period function of global center of system (1.1) is monotonous and
potential energy $G$ is even then the function $x\mapsto\,\frac{G(x)}{x^{2}}$
is monotonous in $(0,+\infty)$.

Parameterizing the periodic orbit by the abscissa $\xi$ of intersection point
between oval $\Gamma$ with positive half $x-$axis, denoted by $\Gamma_{\xi}$.
Concerning the bifurcation of critical periods from the linear isochronous
center, Variation Lemma 3.2 of [10] stated as follows proposes a necessary
condition to the persistence of $2\pi-$periodic solutions.

######  Lemma 2.5.

Let $(x,y)\rightarrow\,(f(x,y),g(x,y))$ be analytic mapping, consider the one
parameter family of perturbed system from the linear center as follows

$\frac{dx}{dt}=-y+\lambda\,f(x,y),\quad\frac{dy}{dt}=x+\lambda\,g(x,y),$

which has a center at the origin for any parameter $\lambda\in\,J$ where real
closed interval $J$ contains $\lambda=0$. Then period function
$T(\xi,\lambda)$ holds for $T_{\lambda}^{\prime}(\xi,0)=0$ if and only if the
integral

$\int_{0}^{2\pi}[xg(x,y)-yf(x,y)]_{x=\xi\cos\theta,y=\xi\sin\theta}\,d\theta=0.$

proof. The conclusion (i) is proved in [10] by using variation of parameter.
Note that the origin $O$ is an linear isochronous center as $\lambda=0$ and
the period function $T(\xi,0)\equiv\,2\pi.$ By the continuous dependence on
the initial value and parameter associated with the solution of differential
equation, the period of closed orbit $\Gamma_{\xi,\lambda}$ is the variation
of $2\pi$, denoted by $T(\xi,\lambda)=2\pi+\rho(\xi,\lambda)$. It is clear
$\rho(\xi,0)=0$ and $\frac{\partial{\rho(\xi,0)}}{\partial{\lambda}}=0$ is
equivalent to vanishing of the following integral

$\int_{0}^{2\pi}[xg(x,y)-yf(x,y)]_{x=\xi\cos\theta,y=\xi\sin\theta}\,d\theta=0.$

That ends the proof.

The integral of perturbation terms
$\lambda[f(x,y)\frac{\partial}{\partial\,x}+g(x,y)\frac{\partial}{\partial\,y}]$
at the cycle passing through $(\xi,0)$ vanishes which proposes a necessary
condition to persistence of $2\pi$ periodic solutions after perturbations.

##  3\. The proof of Theorem 1.1

In this section we intend to prove Theorem 1.1 by twofold: the first one is
the existence of critical periods, another one is the unboundedness of number
of critical periods.

Firstly it is obvious when $\lambda\neq\,0$ the origin of system (1.2) is non-
isochronous center as shown in Lemma 2.3.

Using Lemma 2.1 we get the asymptotic expansion of $T(h,\lambda)$ for $h$
sufficiently small is

(3.1)
$T(h,\lambda)=\frac{2\pi}{\sqrt{1+\lambda}}+\frac{\lambda\,\pi}{4(1+\lambda)^{\frac{5}{2}}}\,h+O(h^{2}).$

It follows from Lemma 2.2 that the system (1.2) is semilinear, hence period
function $T(h)$ for $h\rightarrow\,\infty$ satisfies

$0<\liminf\limits_{h\rightarrow\,+\infty}T(h)\leq\,\limsup\limits_{h\rightarrow\,+\infty}T(h)<+\infty.$

By direct computation we get for system (1.2) that

$\frac{G(x)}{x^{2}}=\frac{1}{2}+\frac{\lambda(1-\cos(x))}{x^{2}},$

its derivative with respect to $x$ is

$[\frac{G(x)}{x^{2}}]^{\prime}(x)=\lambda\frac{x\sin(x)-2+2\cos(x)}{x^{3}}.$

An easy computation shows that function $\frac{G(x)}{x^{2}}$ is monotonous in
$(0,2\pi)$, but not monotonous in $(0,\infty)$, it follows from Lemma 2.4 the
period function of system (1.1) is not monotonous in $(0,+\infty)$. Therefore
potential system (1.2) has critical periodic orbits in the unbounded period
annulus $R^{2}$.

By Lemma 2.5, we intend to consider the variational problem on the period
$T(\xi,\lambda)$. By change of polar coordinates, we get for the
$2\pi-$periodic orbit $\Gamma_{\xi}$ of system (1.2) it satisfies

$\displaystyle\int_{0}^{2\pi}\sin(\xi\cos(s))\cos(s)\,ds$ $\displaystyle=$
$\displaystyle 4\int_{0}^{\frac{\pi}{2}}\sin(\xi\sin(u))\sin(u)\,du$
$\displaystyle=$ $\displaystyle
4\int_{0}^{1}\sin(\xi\,w)\frac{w}{\sqrt{1-w^{2}}}\,dw=2\pi\cdot\,J(1,\xi)=0.$

Where the notation $J(1,\xi)$ stands for the Bessel function of the first kind
[1] which satisfies the following Bessel’s equation of order one

$\xi^{2}\frac{d^{2}\,x}{d\xi^{2}}+\xi\frac{dx}{d\xi}+(\xi^{2}-1)x=0.$

Bessel Function

$J(1,\xi)=\sum\limits_{m=0}^{+\infty}\frac{(-1)^{m}}{m!(m+1)!}(\frac{\xi}{2})^{2m+1}=\frac{\xi}{2}-\frac{\xi^{3}}{16}+\frac{\xi^{5}}{384}-\frac{\xi^{7}}{18432}+\cdots$

has infinitely many simple zeros in $(0,\infty)$ and satisfies $J(1,0)=0$ when
$\xi\rightarrow\,+\infty$, there holds the asymptotical equality

$J(1,\xi)=\sqrt{\frac{2}{\pi\xi}}\cos(\xi-\frac{3\pi}{4})+O(\xi^{-\frac{3}{2}}).$

It implies the period function of center for system (1.2), denoted by
$T(\xi,\lambda)$ which satisfies $h=G(\xi)$, it holds

$T_{\lambda}^{\prime}(\xi,0)=0.$

Note that $T^{\prime}(\xi)=G^{\prime}(\xi)T^{\prime}(h).$ Hence for parameters
$\xi$ and $h$ the monotonicity of the period function $T$ is the same. Thus
the period function $T(\xi,\lambda)$ oscillates and intersects the line
$T(\xi,0)=2\pi$ infinitely many times, thus it has an infinite number of roots
$h_{k}$ such that there holds $T(h_{k},\lambda)=2\pi$ for all
$\lambda\in\,[-1,+\infty)$.

Without lose of generality assume that $\lambda>0$ in the following
discussion, the case $\lambda\in\,[-1,0)$ can be treated similarly.

In the phase plane taking the intersection points $(x_{k},0)=(2k\pi,0)$
between the oval $\Gamma_{h}$ and $x-$axis, the corresponding energy level is
$h_{k}=G(x_{k})=2(k\pi)^{2}$ where $k$ is a positive integer. Due to the
symmetries when $G$ is even, if $\lambda\in[0,+\infty)$ the period of oval
$\Gamma_{h_{k}}$ is given by

(3.2)
$T(h_{k},\lambda)=2\sqrt{2}\int_{0}^{2k\pi}\frac{dx}{\sqrt{h_{k}-G(x)}}\geq\,2\sqrt{2}\int_{0}^{2k\pi}\frac{dx}{\sqrt{h_{k}-\frac{x^{2}}{2}}}=2\pi.$

Similarly consider the motion along periodic orbit passing through points
$(\bar{x}_{k},0)=(2k+1)\pi,0)$, the corresponding energy level is
$\bar{h}_{k}=\frac{(2k+1)^{2}\pi^{2}}{2}+2\lambda$. The least period of the
motion along the period oval is represented as

(3.3) $\displaystyle T(\bar{h}_{k},\lambda)$ $\displaystyle=$ $\displaystyle
2\sqrt{2}\int_{0}^{(2k+1)\pi}\frac{dx}{\sqrt{\bar{h}_{k}-G(x)}}$
$\displaystyle=$ $\displaystyle
2\sqrt{2}\int_{0}^{(2k+1)\pi}\frac{dx}{\sqrt{\frac{(2k+1)^{2}\pi^{2}}{2}-\frac{x^{2}}{2}+\lambda(1+\cos(x))}}$
$\displaystyle\leq$
$\displaystyle\,2\sqrt{2}\int_{0}^{(2k+1)\pi}\frac{dx}{\sqrt{\frac{(2k+1)^{2}\pi^{2}}{2}-\frac{x^{2}}{2}}}=2\pi.$

Combining with expressions (3.2) and (3.3), it can be deduced by Intermediate
Value Theorem that there exists at least one energy level
$h_{i}\in\,(h_{k},\bar{h}_{k})$ such that the period $T(h_{i},\lambda)=2\pi$
for all $\lambda>-1.$ It implies that the conservative system (1.2) has
countable infinitely many $2\pi$ periodic solutions. It follows from Rolle’s
Theorem that completes the proof of Theorem 1.1.

## References

  * [1] G. E. Andrews, R. Askey, R. Roy; Special functions, Encyclopedia of mathematics and its applications (Book 71), Cambridge University Press; 1 edition (January 29, 2001).
  * [2] V. I. Arnold; Loss of stability of autooscillations near a resonance and versal deformations of equivariant vector fields. Funkts Anal. 11 (1977), 1–10.
  * [3] J. Chavarriga, M. Sabatini; A survey of isochronous centers, Qual. Theory Dyn. Syst. 1(1999), 1–70.
  * [4] C. Chicone; Geometric Methods for Two-Point Nonlinear Boundary Value Problems, J. Differential Equations 72 (1988), 360–407.
  * [5] Z. Opial; Sur les périodes des solutions de l’équation différentialle $x^{\prime\prime}+g(x)=0$, Ann. Polon. Math. 10 (1961), 49–72.
  * [6] M. Urabe; Potential forces which yield periodic motions of a fixed period, Jour. Math. Mech. 10 (1961), 569–578.
  * [7] S. N. Chow, J. A. Sanders, On the number of critical points of the period, J. Differential equations 64 (1986) 51–66.
  * [8] C. Chicone, F. Dumortier; Finiteness for critical periods of planar analytic vector fields, Nonlinear Anal. 20 (1993), 315–335.
  * [9] M. Sabatini, Period function’s convexity for Hamiltonian centers with separable variables, ANNALES POLONICI MATHEMATICI 85.2 (2005) 153–163.
  * [10] C. Chicone, M. Jacobs; Bifurcation of critical periods for plane vector fields, Trans. Amer. Math. Soc. 312 (1989), no. 2, 433–486.
  * [11] J. P. Franco̧ise; The successive derivatives of the period function of a plane vector field, J. Differential Equations 146 (1998), 320–335.
  * [12] E. Freire, A. Gasull and A. Guillamon; Period function for perturbed isochronous centres, Qual. Theory Dyn. Syst. 3 (2002), 275–284.
  * [13] A. Gasull, C. Liu, J. Yang; On the number of critical periods for planar polynomial systems of arbitrary degree, J. Differential Equations 249 (2010), no. 3, 684–692.
  * [14] M. Grau and J. Villadelprat; Bifurcation of critical periods from Pleshkan’s isochrones, J. London Math. Soc. 81 (2010), no. 2, 142–160.
  * [15] Ding, T., Huang, H., Zanolin, F.: A priori bounds and periodic solutions for a class of planar systems with applications to Lotka-Volterra Equations, Disc. Cont. Dyn. Syst. 1, 103–117 (1995)
  * [16] L. Yang, X. Zeng; The period function of potential systems of polynomials with real zeros, Bull. Sci. Math. 133 (2009), 555–577.
  * [17] Z. Zhang; A theorem about the differential equation $\ddot{x}+\mu\sin(x)+x=0$ has and only has $n$ limit cycles at the domain $\mid\dot{x}\mid\leq(n+1)\pi$, China. Sci. 10 (1980), 941–948.