\section{Introduction} \label{S:1}
Since the first astrophysical observation of the “dunkle Materie” by Zwicky in 1933, we know that Dark Matter (DM) makes ~27\% of the matter-energy budget in the Universe. Its elemental composition however remains one of the greatest mysteries in all of physics. We recall that DM plays a major role in structure formation in the Universe. There is no lack of DM candidates, and therefore a large number of experimental searches are being undertaken to answer this physics question of fundamental importance. With the direct DM searches being in the forefront, also accelerator experiments (e.g. LHC experiments and other) and also other branches have targeted this question aiming to produce or detect some DM particle. Notably, the ultimate identification of DM species can be only with those who fill the Universe and are also expected in our neighbourhood.
In a recent investigation of the stratospheric temperature  excursions occurring annually around December-January, it was suggested that their planetary dependence points at streams of invisible matter from the dark universe as being the cause \cite{Zioutas:2020PDU}. In this work we go a step further and suggest how to detect them \textit{in situ}, i.e., at the stratosphere where the hours- or day-long heating-ups occur. As it was explained in \cite{Sirks:2020arxiv,Barak:2020arxiv,Zioutas:2020arxiv,Zioutas:2020PDU,Essig:2019arxiv,Meier:1997AP,Henry:2015ApJ}, the widely mentioned axions or WIMPs come not in question, at least not those which various experimental groups are searching for since some decades, since their interaction cross section with normal matter is far too small to account for the observed heating. Therefore, inspired by the reasoning about the origin of the stratospheric temperature anomalies, the stratosphere might be the place with occasionally enhanced converted DM constituents. We outline here a direct experimental search for DM candidates in the stratosphere, i.e., some 40 km above Earth’s surface. This novel idea of direct detection brings the experiment from underground or ground to the upper atmosphere.

In this proposal we focus primarily on dark photons, which are also DM candidates. The reasoning behind this choice is the following: dark photons can convert kinetically to ordinary photons, by contrast to the case with axions or axion-like particles, where a magnetic or electric field is needed to assist their photon conversion [3]. Notably, if the charge density has a plasma frequency which matches the rest mass of the dark photon, the conversion is enhanced due to resonance/coherence effects [3]. Following previous astrophysical observations combined with the axion or axion-like exclusion by the CAST experiment \cite{CAST:2017nat}, the dark photon rest mass might be up to the neV range. The corresponding plasma density is of the order of $10^{4 \pm 2}$ electrons/cm$^{3}$. Interestingly, such electron densities exist above the stratosphere (altitude above 100 km) [2]. Then, the dark photon scenario could fit-in for the heating up. Low speed streams of dark photons incident on the solar system experience solar and planetary gravitational self-focusing towards the Earth, before entering the ionosphere with the appropriate plasma density  some 100s to 1000s km above the Earth’s surface, where they can convert to real photons. The photons from resonant conversion would be RF photons in the MHz range. However, non-resonant conversions at other bands could be happening, albeit not the same enhancement as resonant ones. One such a region of interest is UV range around 6-8 eV because of the absorption edge in the upper stratosphere, which strongly prevents these photons from reaching ground level.

Thus, once an instrument, with its field-of-view pointing away from the Earth, observes temporally enhanced radio signals or UV photons in the stratosphere this will be the new direct DM signature, e.g., for converted dark photons or other ways indirectly. If a detected signature correlates with the local stratospheric temperature, and also exhibits matching planetary dependence, this strongly points to an invisible matter origin of the signal and at the same time excluding background sources.

Once such measurements in the stratosphere have been performed, they also could be correlated with missions in space including the International Space Station.

Such measurements could become the harbinger of a potential major direct discovery of DM constituents. Therefore, a multi-purpose detector system has the potential to unravel more species from the dark sector, which might well show enhanced fluxes in the upper atmosphere.

In short, we propose novel searches for DM constituents, in a location complementary to current DM searches and exploiting new signatures.

The widely assumed annual oscillation of the DM wind is here overruled by planetary dependence \cite{Bertolucci:2017pdu}, which has not been so far considered in direct DM searches.

\section{The concept}
\label{S:2}
We suggest to directly search for streaming DM incident on the Earth’s atmosphere. Our motivation stems also from the fact that a possible signature of DM might be shielded by the atmosphere, in which case ground and underground DM search experiments would be quasi-blind. We propose as the detector location the stratosphere due to the planetary dependence of the temperature already observed there \cite{Zioutas:2020PDU,Essig:2019arxiv} and to the natural occurring plasma in the ionosphere above, which could produce coherent conversion effects.

Furthermore, it is known that solar photon irradiation of the atmosphere shows various characteristic absorption layers. The stratosphere around 40-45\,km above the surface absorbs UV~photons in the energy range $\sim6-8$ \,eV because of the overall oxygen absorption bands \cite{Meier:1997AP}. Thus, any search in this energy band should be done at stratospheric altitudes.

We assume that incoming slow DM constituents first experience gravitational effects by the solar system bodies, including the Sun, and then are transformed into photons in the energy range mentioned above. Only photons emerging from slow DM constituents would show a planetary relationship. On the other hand, light or relativistic cosmic radiation cannot show a planetary dependence, since even the strongest possible gravitational focusing effect by the Sun has a focal length of $\sim 550$ Astronomical Units (A.U.) for photons.

A favored candidate would be the dark photon since it kinetically mixes with real photons without the need for electric or magnetic fields, and, in addition, this conversion could be resonantly enhanced if the plasma density matches the rest mass of the dark photon $(A^{'})$, $\hbar \omega_{plasma} = m_{A'}$ \cite{Gelmini:2020arxiv}.
The interaction of a dark photon with a plasma medium is described for example in \cite{Gelmini:2020arxiv,Dubovsky:2018arxiv}.
The ionospheric plasma density is changing steadily, reaching a peak in electron density at an altitude of about 300 km. The typical electron density of the ionosphere is $N_e=10^{12}$\,electrons/m$^3$ resulting in a plasma frequency of

\begin{equation}
\label{eq:PlasmaFreq}
\omega_{plasma}=2\pi \sqrt{\frac{N_e e^2}{m_e \epsilon_0}}\approx 50 \mbox{[Hz]}\cdot \sqrt{N_e/ \mbox{m}^{-3}}=50\,\mbox{MHz}
\end{equation}

If the dark photon mass is in the neV mass range, corresponding to some MHz in frequency, then the interaction of a dark photon with the ionospheric plasma enhances production of real photons. However, when the dark photon mass is in the eV range or above, the conversion in the ionosphere is off resonance and therefore
smaller \cite{Dubovsky:2018arxiv}.

Radiation flux and power due to dark photon to photon conversion can vary widely depending on kinetic mixing, temporary flux enhancement due to gravitational focusing, resonance conditions, variability of local degree of ionization, and local dark matter density (see \cite{Gelmini:2020arxiv}).

For a $\approx$neV dark photon the coherence effects are determined by level crossings depending on the profile of ionization column as a function of altitude. Using a Chapman electron density profile given by ~\cite{Stankov:2003}

\begin{equation}
\label{eq:elecDensProf}
n_{e}(z )= n_{max} e^{\frac{1}{2}[1-\frac{z-h_{max}}{H}}-e^{(-\frac{z-h_{max}}{H})]}
\end{equation}

with $n_{max}=10^{12} \mbox{electrons} / \mbox{m}^3, H = 100~\mbox{km}, h_{max} = 300~\mbox{km}$, we calculate a scale factor of $|\frac{d\,log( \omega_{plasma}(z)^2)}{dz}|^{-1}\approx 10^4$ for resonant conversion in the the ionosphere~\cite{McDermott:2020arxiv,Garcia:2020}.

For a kinetic mixing value of $\epsilon=10^{-12}$, a dark matter density of $\rho=0.3 \, \mbox{GeV}/\mbox{cm}^{3}$, an antenna effective area of about $1 \, \mbox{m}^2$, and for a very modest gravitational amplification of only 100 from the Sun or solar system bodies, the power received in the MHz frequency range at stratospheric altitude is estimated as $\sim 5\times 10^{-15}$~W~\cite{Horns:2012arxiv}.
The noise at these frequencies is mainly determined by system temperature and bandwidth. The system temperature includes sky noise and receiver noise. For a nominal system temperature of 5000~K, a bandwidth of 1~kHz, the noise power is $P_{noise}\approx 7\times 10^{-17}$~W yielding almost a factor of 100 signal to noise ratio. For weaker signals shielding techniques could be applied to reduce the receiver noise~\cite{Carr:2001}.

For an eV range dark photon the ionospheric plasma is very diluted and the conversions to photons are non-resonant, approximately like those in vacuum. Using power equation from \cite{Gelmini:2020arxiv} we estimate the number of photons from dark photons to be from about $10^{-6}-10^{4}\,\mbox{photons}\cdot\mbox{cm}^{-2}\cdot\mbox{s}^{-1}$.

The sensitivity of an eV range measurement depends on the capabilities of detectors to measure the modulation in the presence of daytime and/or nighttime background.
In the 150-200 nm range the dominant contribution of photon flux is from the Sun. Solar irradiance varies approximately linearly over the mentioned wavelength range from about $0.1\,\mbox{mW}\cdot\mbox{m}^{-2}\cdot\mbox{nm}^{-1}$ to $10 \,\mbox{mW}\cdot\mbox{m}^{-2}\cdot \mbox{nm}^{-1}$ at
the top of the atmosphere, which then, at an altitude of $\approx 45$ km, is reduced by a factor $e$, corresponding to one optical depth ~\cite{Meier:1997AP}.
Different absorption bands in this wavelength region and at these altitudes result in a spectral flux ranging from $10^{8}-10^{11}\,\mbox{photons}\cdot\mbox{cm}^{-2}\cdot\mbox{s}^{-1} \cdot\mbox{nm}^{-1}$~\cite{Meier:1997AP}. For daytime measurements, in order to suppress the solar light, a coronagraph could be used with a typical suppression of $10^6$ resulting in a photon flux of $5\times 10^{3} -5\times 10^{6} \,\mbox{photons}\cdot\mbox{cm}^{-2}\cdot\mbox{s}^{-1}$ over this range of wavelengths.  If daytime measurement are made, this is the main photon background.

Missions have made night-time measurements in low Earth orbit ($\approx 600 \, \mbox{km}$) in this energy range \cite{Henry:2015ApJ}.  They report a baseline spectral flux of $\approx 10^{4}\,\mbox{photons}\cdot\mbox{cm}^{-2}\cdot\mbox{s}^{-1}\cdot\mbox{nm}^{-1}$ which is primarily from diffuse galactic and extragalactic light and airglow of the atmosphere.  Reducing the flux by a factor of one optical depth to scale it at altitudes of 45 km, integrating over the solid angle, and considering different absorption bands in this wavelength region, the estimated cosmic noise level is $500-3\times 10^{5}\,\mbox{photons}\cdot\mbox{cm}^{-2}\cdot\mbox{s}^{-1}$ for the 150 nm -200 nm wavelength range. This would be the main background for night-time measurements.

The overall signature would be a modulated 6-8\,eV photon flux or an RF power in the MHz range, concurrent with some planetary correlation or simultaneous temperature excursions in the stratosphere. Thus, at the beginning, the primary detector choices would be UV photon detectors and low noise radio receivers. Other detectors could be added on to probe other energy bands, wider fields of view, direction discrimination, and ancillary parasitic detectors as payload limit allows.

Small payloads could be sent out first in weather balloon flights lasting hours to measure the day and night background and establish the best detector combinations for these measurements. This could then be followed by bigger and longer duration balloon missions with multi-purpose detectors.

\begin{figure}[h]
\centering\includegraphics[width=1.0\linewidth]{DarkPhotonStream_v3.jpg}
\caption{Cartoon illustration of the direct search concept (not to scale). A gravitationally focused stream of dark photons partially converts into real photons in a region farther that $\approx$100\,km above the surface of the earth, where the plasma density has the resonant value. Converted photons are then absorbed in the upper stratosphere, $\approx$40\,km above the surface, causing the observed local temperature increase \cite{Zioutas:2020PDU}. A photon detector placed in the upper stratosphere could directly measure excess photons coming from converted dark photons (A').}
\end{figure}

\section{Possible detectors}
\label{S:3}

Balloon flights carrying a payload of several hundred kilograms in the stratosphere can last for days, providing the opportunity to attempt direct detection of DM as outlined above.
The ideal device fitting the proposed search would be a highly efficient, wide sensitive area, low noise photon detector covering a 4$\pi$ Field-Of-View (FOV), divided at least in two separate halves. The detector sections pointing away from Earth would search for variations in the incoming photon flux in the energy range of interest. The opposite sectors would act as background monitors.
For the detection of 6-8 eV photons from converted dark photons, a promising option would be using Silicon Photomultipliers (SiPMs). These detectors are solid-state devices,
Commercially available SiPMs can be arranged in arrays to cover active areas of several cm$^2$, and identical arrays could be combined with light guides to effectively cover large opposite portions of the sky. This would allow the pairwise comparison between detectors looking skywards and earthwards.
A preliminary test payload could also be equipped with UV sensitive photo-diodes for complementary monitoring.
In a possible payload configuration, SiPMs could be combined with RF receivers, such as those developed for radio astronomy, to provide sensitivity in the MHz range. For instance receivers with a dipole-like antenna developed for The Long Wavelength Array~\cite{Henning:2010arxiv} would suffice.
Balloon payloads could also be complemented with unconventional state-of-the-art detectors such as ultra-high sensitivity opto-mechanical force sensors \cite{Karuza:2016pdu}. The force sensors could be operated either in the dark mode where it would look for interactions with dark sector particles with its surface, or light mode where also the radiation pressure from photons will be sensed \cite{Arguedas:2019pdu}.

\section{Conclusions}
\label{S:4}

We have presented here a novel proposal to search for Dark Sector constituents converting into Standard Model particles, which then cause observable effects at upper stratospheric altitudes. The proposal relies on the idea of placing detectors in the stratosphere.

Although we focus on the example of dark photons kinetically converting into real photons, where the mixing can also be resonantly enhanced by the intervening ionospheric plasma, the search could be extended to other candidates from the Dark Sector by envisioning a suitable multi-purpose detector. The proposed search could also include instruments in orbit.

The other important feature is that the signature of Dark Sector conversion must show a dependence on the relative positions of the planets in our solar system. Conventional cosmic rays are mostly relativistic and therefore the solar system cannot exert a significant gravitational impact on their trajectories, certainly not on photons.
The new search strategy proposed in this work is fully complementary to underground experiments. Even surface detectors could not perform this type of search for some Dark Sector candidates as they may be shielded by the Earth’s atmosphere with about $10^3 \, \mbox{g} \cdot \mbox{cm}^{-2}$.

\section{Acknowledgements}
\label{S:5}
For M. Maroudas, this research is co-financed by Greece and the European Union (European Social Fund- ESF) through the Operational Programme ‘‘Human Resources Development, Education and Lifelong Learning’’ in the context of the project ‘‘Strengthening Human Resources Research Potential via Doctorate Research’’ (MIS-5000432), implemented by the State Scholarships Foundation (IKY). This work has been partially supported by University of Rijeka, Croatia, grant number 18-126.

\bibliographystyle{elsarticle-num-names}


\section{Introduction\label{secI}}
{The teleparallel equivalent of general relativity (TEGR) is an equivalent formulation of Einstein's general relativity (GR), in that, the formulation ensures equivalence of the field equations as well as of the test particle equations of motion. For a review see the paper by Maluf~\cite{Maluf} where an account of the history of teleparallel theories of gravity is given. The two theories, TEGR and GR, use different connections, the curvature-less Weitzenb\"{o}ck connection~\cite{W} and the torsion-less Levi-Civita connection~\cite{LC}, respectively. The Weitzenb\"{o}ck connection has the property that it allows for a definition of a condition for absolute parallelism in space–time~\cite{Maluf}, hence the given name of ``teleparallel''. The tensors and invariants associated to TEGR exhibit no curvature but torsion only, that is, the information concerning the  gravitational field effects are encoded in the torsion tensor instead of the Riemann tensor, as it is the case in GR.}

{Despite their equivalence, the two theories are conceptually different. All physical features and results known in GR are also
described in the TEGR. The converse is not true: The TEGR approach allows the consideration of additional concepts and definitions. In the TEGR one meets tensors with three indices some of which are behind the definition of the gravitational energy–momentum which is consistent with the field equations~\cite{19,20}. The concept and definition of the gravitational angular momentum are also introduced in the TEGR~\cite{Maluf}. From the point of view of test particles, there is no notion of geodesics in the TEGR but only force equations and the source of force is torsion, more precisely the force is sum of a torsion-tensor component times two components of the velocity vector~\cite{book}. In the TEGR it is also possible to introduce the concept of inertia-free frames~\cite{book} and to split particle dynamical effects into distinct gravitational effects and purely inertial effects, that is, it is possible to separate inertia from gravity~\cite{nonzero}. Another example of conceptual difference between the two theories is provided in this work following Eq.~\eqref{pr}.}

{We have already mentioned some motivations for the investigation of the TEGR: Introduction of new concepts and unambiguous definitions such as that of the gravitational energy–momentum and other definitions~\cite{Maluf} and the clear distinction of purely inertial effects from gravitational ones. Other motivations for the TEGR is that the theory admits some extensions by adding quadratic and higher-order torsion terms to the action making it a good cosmological dark energy model without truly adding exotic matter to the cosmological field equations. Said otherwise, these extra terms added to the action, the so-called $f(T)$ gravity, have their counterparts in the cosmological field equations playing the ``role'' of exotic matter~\cite{Maluf,Maluf2}. The extended TEGR may provide a theoretical model to the late time universe acceleration problem~\cite{4}.}

Some rotating cylinders in general relativity sourced by anisotropic fluids and $\Lambda =0$ have been determined in~\cite{cz,lb,w4}. The purpose of this work is to consider the theories of GR and TEGR coupled to a cosmological constant and determine cylindrically symmetric rotating wormholes and other solutions sourced by anisotropic fluids.

Determining rotating and static solutions around an infinite axis is still a reviving topic. In GR there is a set of rotating cylindrically symmetric perfect fluid solutions which may be appropriate as matched interiors. Among the known solutions in GR we find the rotating dust of Vishveshwara and Winicour~\cite{fd1}, the perfect fluid sources with non-zero pressure of da Silva et al.~\cite{Lewis2}, Davidson~\cite{fd3,Davidson} and Ivanov~\cite{fd5}, and the family of Krasi\'{n}ski~\cite{fd6}. The solutions that will be constructed in this work share some physical and geometrical properties with the solutions known in the literature and have some other new properties. Among the new properties we mention that their angular velocity is proportional to the cosmological constant.

Investigating wormhole solutions is another, {rather many-fold}, reviving topic. {First of all, wormholes are special types of solutions to the field equations of gravity theories which contain two distant asymptotic regions or sheets, with a throat connecting the two and providing a shortcut for long journeys from one asymptotic region to another~\cite{MT,V}. Their observation has not been confirmed yet but they may be very well lurking in the universe. In Ref.~\cite{types} it was shown that the diameter of the shadow of type I supermassive wormhole overlaps with that of the black hole candidate at the center of the Milky Way. This shows that the existing up-to-date millimeter-wavelength very long baseline interferometry facilities do not lead to differentiate the supermassive black hole candidate at the center of the Milky Way from a possible type I supermassive wormhole.} Very recently, it was shown that the active galactic nuclei exhibit wormhole behaviors rather than supermassive black holes due to their gamma radiation resulting from collision of accreting flows~\cite{agn}. {Said otherwise, wormholes may be spotted in the sky upon detecting the gigantesque display of gamma rays that results from the collision of matter coming out of one mouth of the wormhole with infalling matter~\cite{agn}. A nearly similar conclusion was made in Ref.~\cite{types}: ``Other signals from the galaxy, as the motion of orbiting hot spots, may lead to draw a conclusion concerning the nature of the candidate''}.

{The interest to wormholes goes back to 1935. Einstein and N. Rosen were the first who described how two distant regions of spacetime can be joined together, creating a bridge between them. Most of the known solutions, and their list is too large to be cited here, are endowed with spherical symmetry, until recently a couple of wormhole solutions endowed with cylindrical symmetry were determined (see~\cite{symmetry} are references therein). The construction of analytical wormhole solutions, with spherical or cylindrical symmetry, provides the scientific community with theoretical tools for further investigations (calculations of shadow, quasi-normal modes, quasi-periodic oscillations, etc), as was done in~\cite{types,hot1,hot2}, and for computer simulations, as was done in~\cite{col1,col2} to study the collisional processes in the geometry of rotating wormholes.}

In Sec.~\ref{sectegr} we introduce the mathematical tools needed for the TEGR along with the necessary field equations. In Sec.~\ref{secT} we reduce the field equations. In Sec.~\ref{secEOS} we restrict ourselves to anisotropic fluids with $p_r=\omega_r \rho$, $p_\phi=\omega_\phi \rho$ and $p_z=\omega_z \rho$, where the equation-of-state (EoS) parameters ($\omega_r,\,\omega_\phi,\,\omega_z$) are constants constrained by $-1\leq \omega_r\leq 1$, $-1\leq \omega_\phi\leq 1$ and $-1\leq \omega_z\leq 1$. Section~\ref{secomega} is devoted to the construction of rotating wormholes and Sec.~\ref{rws} is devoted to the discussion of their physical and geometrical properties. In Sec.~\ref{sectas} we provide cosmological rotating solutions with two axes of symmetry at finite proper distance. Our final conclusions are given in Sec.~\ref{secc}. {An appendix section has been added to provide some useful formulas pertaining to Sec.~\ref{secomega}}.

{{\section{The teleparallel equivalent of general relativity\label{sectegr}}}
In the TEGR~\cite{HS9} the vielbein vector fields $e_a=e_{a}{^{\mu}}\partial_\mu$ are taken as fundamental variables instead of the metric $g_{\mu\nu}$, related to each other by
\begin{equation}\label{i1}
	g_{\mu \nu} =
	\eta_{ab}e^{a}{_{\mu}}e^{b}{_{\nu}},\quad g^{\mu \nu} =
	\eta^{ab}e_{a}{^{\mu}}e_{b}{^{\nu}},\quad e=\sqrt{|g|},
\end{equation}
with $\eta_{ab}=\text{diag}(+1,\,-1,\,-1,\,-1)$ being the metric of the 4-dimensional Minkowski spacetime and $e\equiv |\det(e^{a}{_{\mu}})|=\sqrt{|g|}$. In this work the tetrad indices, $a,\,b,\,c\cdots,k,\,l$, and Greek coordinate indices run from 0 to 3.

The field equation in teleparallel gravity in the presence of matter fields take the form~\cite{Bengochea,FF7}
\begin{align}\label{i2}
& I_\mu{^\nu}:=-\delta^\nu_\mu~\frac{f}{2} +2\Big[e^{-1}{e^a}_\mu\partial_\rho\Big(e{e_a}^\alpha
{S_\alpha}^{\rho \nu}\Big)-{T^\alpha}_{\lambda \mu}{S_\alpha}^{\nu \lambda}\Big]\nonumber\\
& =-\kappa T_{\text{(mat)}\,\mu}{}^\nu,
\end{align}
where $\kappa$ is the gravitational constant,
\begin{equation}\label{i3}
f(T)=T+2\Lambda ,
\end{equation}
and $T_{\text{(mat)}\,\mu}{}^\nu$ is the matter stress-energy tensor (SET) which we assume to be that of an anisotropic fluid of the form
\begin{align}\label{i4}
T_{\text{(mat)}\,\mu}{}^\nu & =\rho u_\mu u^\nu+\sum_{a=1}^3p_a{e_a}_\mu {e_a}^\nu\nonumber\\
& =(\rho+p_1)u_\mu u^\nu-p_1\delta^\nu_\mu+(p_2-p_1){e_2}_\mu {e_2}^\nu\nonumber\\
& +(p_3-p_1){e_3}_\mu {e_3}^\nu .
\end{align}
Here we have chosen ${e_0}^\nu=u^\nu$ to be the four-velocity vector of the fluid. The remaining quantities used in~\eqref{i2}, including the torsion $T$, are defined by\footnote{$S^{\alpha\mu\nu}$ may be given in a more compact form as:$$S^{\alpha\mu\nu}=\frac{1}{4}(T^{\nu\mu\alpha}+T^{\alpha\mu\alpha}-T^{\mu\nu\alpha})-\frac{1}{2}g^{\alpha\nu}{T^{\sigma\mu}}_{\sigma}+\frac{1}{2}g^{\alpha\mu}{T^{\sigma\nu}}_{\sigma}.$$}
\begin{align}\label{i5}
&{T^\alpha}_{\mu \nu}={e_b}^\alpha
(\partial_\mu{e^b}_\nu-\partial_\nu{e^b}_\mu),\nonumber\\
&K_{\alpha\mu\nu}=\frac{1}{2}~(T_{\mu\alpha\nu}+T_{\nu\alpha\mu}-T_{\alpha\mu\nu}),\nonumber\\
&S^{\alpha\mu\nu}=\frac{1}{2}~(K^{\mu\nu\alpha}-g^{\alpha\nu}{T^{\sigma\mu}}_{\sigma}+g^{\alpha\mu}{T^{\sigma\nu}}_{\sigma}),\nonumber\\
&T=T_{\alpha\mu\nu}S^{\alpha\mu\nu}.
\end{align}
In the definition of ${T^\alpha}_{\mu \nu}$ the connection $\omega^a{}_{b\mu}$ has been set equal to zero as this is always possible in teleparallel gravity~\cite{cov}.

In cylindrical coordinates ($x^0=t$, $x^1=r$, $x^2=\phi$, $x^3=z$), we introduce the following non-diagonal vielbein to describe rotating solutions
\begin{equation}\label{i6}
\big(e^{a}{_{\mu}}\big)= \begin{pmatrix}
{\rm e}^{\gamma (r)} & 0 & -{\rm e}^{-\gamma (r)} \Omega (r) & 0 \\
0 & {\rm e}^{\alpha (r)} & 0 & 0 \\
0 & 0 & {\rm e}^{\beta (r)} & 0 \\
0 & 0 & 0 & {\rm e}^{\mu (r)}
\end{pmatrix},
\end{equation}
resulting in the metric
\begin{multline}\label{m}
{\rm d}s^2=  {\rm e}^{2\gamma(r)}[ {\rm d}t - \Omega(r){\rm e}^{-2\gamma(r)}{\rm d}\phi ]^2- {\rm e}^{2\alpha(r)}{\rm d}r^2 \\- {\rm e}^{2\beta(r)}{\rm d}\phi^2
- {\rm e}^{2\mu(r)}{\rm d}z^2.\qquad\qquad\qquad\qquad
\end{multline}
Using all that in~\eqref{i5} we obtain
\begin{multline}\label{T}
T=\frac{1}{2} {\rm e}^{-2 (\alpha + \beta + \gamma )} (\Omega '-2 \Omega  \gamma ')^2\\
+2 {\rm e}^{-2 \alpha } (\beta ' \gamma '+\beta ' \mu '+\gamma ' \mu ').
\end{multline}
The nonvanishing components of the SET~\eqref{i4} are
\begin{align}\label{i4nv}
&T_{\text{(mat)}\,t}{}^t=\rho,\quad T_{\text{(mat)}\,r}{}^r=-p_r,\quad T_{\text{(mat)}\,\phi}{}^\phi=-p_\phi,\nonumber\\
&T_{\text{(mat)}\,z}{}^z=-p_z,\quad T_{\text{(mat)}\,\phi}{}^t=-\Omega (\rho+p_\phi){\rm e}^{-2 \gamma },
\end{align}
where we have set $p_1=p_r,\,p_2=p_\phi,\,p_3=p_z$.

It has become customary to introduce the vortex $\omega(r)$, which is the norm of the curl of the tetrad $e_{a}{^{\mu}}$~\cite{Bonnor} (see also~\cite{w1,w2,w3,w4}). This is related to $\Omega(r)$ by\footnote{The curl $\omega^\mu$ of $e_{a}{^{\mu}}$ is \[\omega^\mu:=\frac{1}{8}~\epsilon^{\mu\nu\rho\sigma}e_{a\nu}e^{a}{_{\rho;\sigma}},\] yielding $\omega^t=\omega^r=\omega^\phi=0$, $\omega^z=-{\rm e}^{\gamma-\alpha-\beta-\mu}(\Omega {\rm e}^{-2\gamma})'/2$ and $\omega=\sqrt{\omega^\mu\omega_\mu}={\rm e}^{\gamma-\alpha-\beta}(\Omega {\rm e}^{-2\gamma})'/2$.}
\begin{equation}\label{w}
\Omega (r):=2 {\rm e}^{2 \gamma (r)} \int^r {\rm e}^{\alpha (x)+\beta (x)-\gamma (x)} \omega (x){\rm d}x ,
\end{equation}
yielding
\begin{equation}\label{Ts}
T=2 \omega ^2+2 {\rm e}^{-2 \alpha } (\beta ' \gamma '+\beta ' \mu '+\gamma ' \mu ').
\end{equation}

One may bring the field equations~\eqref{i2} to\footnote{In this work {\[R^{\alpha}{}_{\beta\mu\nu}:=\partial_\nu\Gamma^{\alpha}_{\beta\mu}-\partial_\mu\Gamma^{\alpha}_{\beta\nu}+\Gamma^{\alpha}_{\eta\nu}\Gamma^{\eta}_{\beta\mu}-\Gamma^{\alpha}_{\eta\mu}\Gamma^{\eta}_{\beta\nu}.\]}}
\begin{equation}\label{eff}
G_{\mu}{}^\nu=-\kappa \tau_{\mu}{}^\nu
\end{equation}
where $G_{\mu}{}^\nu$ is the Einstein tensor and $\tau_{\mu}{}^\nu$ is the total SET including the cosmological constant and is defined by
\begin{align}\label{tau}
	&\tau_{t}{}^t=\rho-\frac{\Lambda}{\kappa },\quad \tau_{r}{}^r=-p_r-\frac{\Lambda}{\kappa },\nonumber\\
	&\tau_{\phi}{}^\phi=-p_{\phi }-\frac{\Lambda}{\kappa },\quad \tau_{z}{}^z=-p_z-\frac{\Lambda}{\kappa },\nonumber\\
	&\tau_{\phi}{}^t=-\Omega  (\rho +p_{\phi }){\rm e}^{-2 \gamma }.
\end{align}
Since $\nabla_\nu G_{\mu}{}^\nu\equiv 0$, one must have
\begin{equation}\label{must}
\nabla_\nu \tau_{\mu}{}^\nu= 0.
\end{equation}
The solutions we will derive in this work will satisfy the field equations of GR~\eqref{eff} and of TEGR~\eqref{i2}. For TEGR, the solutions are constrained by~\eqref{Ts}.

\section{Reducing the field equations\label{secT}}
Given that $T_{\text{(mat)}\,t}{}^\phi=0$ and $\tau_{t}{}^\phi=0$, the line $I_t{^\phi}=0$ in~\eqref{i2}, and the component  $_{t}{}^\phi$ of Eq.~\eqref{eff}, reduce upon using~\eqref{i4nv} and~\eqref{tau} to \[I_t{^\phi}=G_{t}{}^\phi=-{\rm e}^{\gamma- \alpha -\beta}[\omega(2\gamma'+\mu')+\omega']=0,\] yielding
\begin{equation}\label{omega}
\omega(r)=\omega_0{\rm e}^{-2\gamma -\mu},
\end{equation}
where $\omega_0$ is a constant of integration. Using this in~\eqref{w} we obtain
\begin{equation}\label{w2}
\Omega (r)=2\omega_0 {\rm e}^{2 \gamma (r)} \int^r {\rm e}^{\alpha (x)+\beta (x)-3\gamma (x)-\mu(x)} {\rm d}x .
\end{equation}

The expression of $G_{r}{}^r$ is just half that of $T$~\eqref{Ts}, $G_{r}{}^r=\omega ^2+{\rm e}^{-2 \alpha } (\beta ' \gamma '+\beta ' \mu '+\gamma ' \mu ')=T/2$, and this implies using~\eqref{eff} and~\eqref{tau} that $p_r$ is given by
\begin{equation}\label{pr}
p_r= \frac{T}{2 \kappa}-\frac{\Lambda}{\kappa} ,
\end{equation}
where the last term proportional to $\Lambda$ is the radial pressure due to the cosmological constant and the first term, $T/(2\kappa)$, is the radial pressure generated by a constant torsion. Thus, in the TEGR, a {\textquoteleft nonvanishing\textquoteright} torsion generates a nonvanishing radial pressure $T/(2\kappa)$ while in GR the relation~\eqref{pr} is written as $G_{r}{}^r=\kappa p_r+\Lambda$ and it merely expresses the fact that the geometric entity $G_{r}{}^r$ is proportional to the sum of the pressures $p_r$ and $\Lambda/\kappa$. {This is another example of conceptual difference between the two theories that we mentioned in the Introduction. This is at the level of the field equations, where the torsion appears as a force while its counterpart in GR, the curvature scalar $\mathcal{R}$, does not assume a dynamical role in the field equations.}

Substituting~\eqref{tau} into~\eqref{must} we obtain
\begin{equation}\label{ntau}
-4 p_r'-(\beta '+\gamma '+\mu ') p_r-\gamma ' \rho +\beta ' p_{\phi }+\mu ' p_z=0 .
\end{equation}

\section{Anisotropic fluids\label{secEOS}}
We consider the simple case where  $p_r=\omega_r \rho$, $p_\phi=\omega_\phi \rho$ and $p_z=\omega_z \rho$, with the EoS parameters ($\omega_r,\,\omega_\phi,\,\omega_z$) being constants generally constrained by
\begin{equation}\label{obey}
-1\leq \omega_r\leq 1,\quad -1\leq \omega_\phi\leq 1,\quad -1\leq \omega_z\leq 1 .
\end{equation}
The differential equation~\eqref{ntau} becomes
\begin{equation}\label{rho}
(\omega _r-\omega _{\phi }) \beta '+(1+\omega _r) \gamma '+(\omega _r-\omega _z) \mu '+4 \omega _r \frac{\rho
	'}{\rho }=0,
\end{equation}
resulting in
\begin{equation}\label{wr}
\rho =\rho _0\exp\Big[\frac{(\omega _{\phi }-\omega _r) \beta -(1+\omega _r) \gamma +(\omega _z-\omega _r) \mu}{4 \omega _r}\Big],
\end{equation}
where $\rho_0$ is a constant of integration. Using this in~\eqref{pr} we can evaluate the torsion from
\begin{equation}\label{torsion}
	T=2(\kappa\omega_{r}\rho+\Lambda).
\end{equation}

Since the radial coordinate $r$ can be changed at will by a coordinate transformation $r\to \bar{r}$, from now on we fix the coordinate gauge to be
\begin{equation}\label{gauge}
\alpha =\beta +\gamma +\mu .
\end{equation}
In this gauge the independent field equations emanating from~\eqref{eff}
\begin{equation}\label{R}
R_{\mu}{}^\nu=-\kappa [\tau_{\mu}{}^\nu -(1/2)\delta_{\mu}{}^\nu \tau_{\sigma}{}^\sigma],
\end{equation}
reduce to
\begin{align}
\label{eqs1}&2{\rm e}^{-2 \alpha } \gamma ''+4 \omega ^2=2\Lambda + \kappa  (1+\omega _r+\omega _{\phi }+\omega _z) \rho, \\
\label{eqs2}&2{\rm e}^{-2 \alpha } \mu ''=2\Lambda + \kappa  (-1+\omega _r+\omega _{\phi }-\omega _z) \rho, \\
\label{eqs3}&2{\rm e}^{-2 \alpha } \beta ''-4 \omega ^2=2\Lambda + \kappa  (-1+\omega _r-\omega _{\phi }+\omega _z) \rho,
\end{align}
which are the $_{t}{}^t$, $_{z}{}^z$ and $_{\phi}{}^\phi$  equations~\eqref{R}, respectively. The $_{\phi}{}^t$ equation~\eqref{R} is a combination of~\eqref{eqs1} and~\eqref{eqs3} and the $_{r}{}^r$ equation~\eqref{R} is a combination of~\eqref{eqs1}, \eqref{eqs2}, \eqref{eqs3} and $G_{r}{}^r=\omega ^2+{\rm e}^{-2 \alpha } (\beta ' \gamma '+\beta ' \mu '+\gamma ' \mu ')$ [recall that $G_{r}{}^r=T/2$].

\section{Cosmological rotating wormholes\label{secomega}}
From now on, we assume $\omega_{r}\neq 0$ to ensure that~\eqref{wr} remains valid. We first construct cosmological rotating cylindrical wormholes to GR. The counterpart solutions to TEGR are the same with the torsion given by~\eqref{torsion}.

\begin{figure}[h]
\centering
  \includegraphics[width=0.49\textwidth]{Fig1.pdf} \\
  \caption{\footnotesize{The embedding diagram of the wormhole solution~\eqref{metric2} for $\delta=4/5$, $q=3/5$ and $u_0=1/10$.}}\label{Fig1}
\end{figure}
\begin{figure}[h]
	\centering
	\includegraphics[width=0.49\textwidth]{Fig2a.pdf} \includegraphics[width=0.49\textwidth]{Fig2b.pdf} \includegraphics[width=0.49\textwidth]{Fig2c.pdf} \\
	\caption{\footnotesize{The embedding diagram of the wormhole solution~\eqref{metric2} showing the upper sheet only for $u_0=1/10$. (a) $\delta=4/5$ and $q=3/5$ (upper panel), (b) $\delta=4/50$ and $q=3/5$ (intermediate panel), (c) $\delta=4/50$ and $q=3/50$ (lower panel).}}\label{Fig2}
\end{figure}

We look for solutions with a constant shift function, that is, $\gamma=0$ and $\mu=0$. In this case $\alpha=\beta$~\eqref{gauge} and Eq.~\eqref{eqs2} implies that $\rho$ is a constant given by $\rho_0=2\Lambda/[\kappa (1-\omega _r-\omega_{\phi }+\omega _z)]$. Since $\rho$ is constant, Eq.~\eqref{wr} implies
\begin{equation}\label{equal}
	\omega _{\phi }=\omega _r,
\end{equation}
and finally
\begin{equation}\label{rho1}
	\rho_0=\frac{2\Lambda}{\kappa (1-2\omega _r+\omega _z)}.
\end{equation}
Equations~\eqref{omega} and~\eqref{eqs1} imply that $\omega^2$ is also a constant proportional to the cosmological constant. This is given by
\begin{equation}\label{o02}
	\omega_0^2=\frac{\Lambda (1+\omega_{z})}{1-2\omega _r+\omega _z}=\frac{\kappa (1+\omega_{z})}{2}~\rho_0.
\end{equation}
Using these vales of ($	\rho_0,\,\omega_0^2$) in~\eqref{eqs3} we bring it to the form (recall $\alpha=\beta$)
\begin{equation}\label{form}
	{\rm e}^{-2 \alpha } \alpha ''=q^2\equiv \frac{2\Lambda (1-\omega _r+2\omega _z)}{1-2\omega _r+\omega _z}.
\end{equation}
This equation can be integrated in all three cases $q^2<0$, $q^2=0$ and $q^2>0$. In the case $q^2<0$ we obtain cosmological rotating solutions with two axes of symmetry at finite proper distance (see Sec.~\ref{sectas}). The case $q^2>0$ (we may assume $q>0$) yields cosmological rotating wormholes where the solution ${\rm e}^{2 \alpha }$  is brought to the form
\begin{equation}\label{alpha}
{\rm e}^{2 \alpha }={\rm e}^{2 \beta }=\frac{c^2\sec^2(c r)}{q^2},
\end{equation}
with $c>0$ being a constant of integration. As to $\Omega$ is obtained from~\eqref{w2}
\begin{equation}\label{O}
	\Omega (r)=2\omega_0\int^r {\rm e}^{2\alpha (x)} {\rm d}x =\frac{2c\omega_0}{q^2}~\tan (cr),
\end{equation}
where we have dropped an additive constant of integration. The metric takes the form
\begin{multline}\label{metric}
	{\rm d}s^2=  \Big[ {\rm d}t -\frac{2c\omega_0\tan (cr)}{q^2}~{\rm d}\phi \Big]^2- \frac{c^2\sec^2(c r)}{q^2}~{\rm d}r^2 \\- \frac{c^2\sec^2(c r)}{q^2}~{\rm d}\phi^2
	- {\rm d}z^2.\qquad\qquad\qquad\qquad
\end{multline}
We see that the spherical radius ${\rm e}^{\beta }=c\sec(c r)/q$ has a minimum value at $r=0$ and increases as $|r|$ increases. On introducing the new radial coordinate $u=c\tan(cr)/q$ and the new constant $u_0^2=c^2/q^2$, we bring the metric~\eqref{metric} to the manifestly wormhole form
\begin{multline}\label{metric2}
	{\rm d}s^2=  \Big( {\rm d}t -\frac{2\omega_0 u}{q}~{\rm d}\phi \Big)^2- \frac{1}{q^2(u^2+u_0^2)}~{\rm d}u^2 \\- (u^2+u_0^2){\rm d}\phi^2
	- {\rm d}z^2.\qquad\qquad\qquad\qquad
\end{multline}
A very similar solution describing a rotating wormhole, which is a solution to the Einstein–Maxwell equations, was determined in~\cite{symmetry}.

The TEGR wormhole solution is also given by~\eqref{metric} and~\eqref{metric2} with $T=2\omega_0^2>0$~\eqref{torsion}, which is a positive constant.

The metric component
\begin{equation}\label{gff}
g_{\phi\phi}=-\Big(\frac{q^2-4\omega_0^2}{q^2}~u^2+u_0^2\Big),
\end{equation}
is manifestly negative if $q^2-4\omega_0^2=2\Lambda(1+\omega_r)/(-1+2\omega_r-\omega_z)\geq 0$ signaling the absence of closed timelike curves (CTCs)~\cite{exact}. For  $q^2-4\omega_0^2< 0$, CTCs occur at large values of $|u|$.

Let us see under which conditions the constraints $q^2>0$, $\omega_0^2\geq 0$ and $\rho_0>0$ are satisfied simultaneously. The gravitational constant $\kappa$ being positive, we obtain the following constraints on  the EoS parameters ($\omega_r,\,\omega_\phi,\,\omega_z$) assuming they obey the general inequalities~\eqref{obey}. If $\Lambda<0$ we obtain
\begin{equation}\label{c1}
\frac{1}{3}<\omega_r\leq 1\;\text{ and }\;\frac{\omega_r-1}{2}<\omega_z<2\omega_r-1 .
\end{equation}
If $\Lambda>0$ we obtain
\begin{align}
\label{c2}&-1\leq \omega_r\leq \frac{1}{3}\;\text{ and }\;\frac{\omega_r-1}{2}<\omega_z\leq 1,\;\text{ or}\\
\label{c3}&\frac{1}{3}<\omega_r< 1\;\text{ and }\;2\omega_r-1<\omega_z\leq 1.
\end{align}
If, however, we want to avoid the presence of CTCs, we have to impose the fourth constraint $q^2-4\omega_0^2\geq 0$. This results in
\begin{equation}\label{c4}
\Lambda>0,\;\;\omega_r=-1\;\text{ and }\;-1< \omega_z\leq 1.
\end{equation}
This shows that, for a positive cosmological constant, it is always possible to have a rotating wormhole with a positive energy density and no CTCs. The anisotropic fluid is isotropic in a plane perpendicular to the axis of rotation with $\omega_{\phi }=\omega_r=-1$.

{To construct the embedding diagram of the metric~\eqref{metric2} we introduce the radial coordinate $R$ defined by $R^2=\delta u^2+u_0^2$ with $\delta\equiv (q^2-4\omega_0^2)/q^2$. For a constant time-slice and $z=\text{const.}$, the two-dimensional spatial metric takes the form \[	{\rm d}s_2^2=\frac{R^2{\rm d}R^2}{q^2(R^2-u_0^2)[R^2-(1-\delta)u_0^2]}+R^2	{\rm d}\phi^2.\] This is to be embedded in the three-dimensional
Euclidean space \[{\rm d}s_3^2={\rm d}R^2++R^2	{\rm d}\phi^2+{\rm d}Z^2.\] The embedded diagram is a surface of revolution, symmetric with respect to the plane $Z=0$, where $Z\equiv Z(R)$ with \[Z(R)=\pm\int_{u_0\sqrt{1+\delta}}^R\sqrt{\frac{x^2}{q^2(x^2-u_0^2)[x^2-(1-\delta)u_0^2]}-1}~{\rm d}x.\] Plots of the surface of revolution as depicted in Figs.~\ref{Fig1} and~\ref{Fig2}. In Fig.~\ref{Fig1} we have shown the two sheets of the wormhole for some set of the parameters ($\delta,\,q,\,u_0$) while in Fig.~\ref{Fig2} only the upper sheet is depicted for different values of the parameters ($\delta,\,q,\,u_0$).}

\section{Physical and geometrical properties of the rotating wormhole solution\label{rws}}
Since the wormhole solutions~\eqref{metric} and~\eqref{metric2} are not asymptotically flat, they are appropriate as matched interiors. Another reason why they are so is that the speed of rotation, $\Omega =2\omega_0 u/q$~\eqref{O}, increases linearly as one moves away from the axis of rotation at $u=0$. Consequently the linear speed of each fluid particle increases unceasingly as one moves away from the axis of rotation. Following the work done in~\cite{symmetry} one can match these wormholes to external rotating flat Minkowskian metrics in such a way that the radii of the cylindrical junction surfaces are chosen so that the interior wormhole region is exempt from CTCs.

Thus, the constraints~\eqref{c1}-\eqref{c3}, for $\Lambda<0$ and $\Lambda>0$, are fairly enough to obtain physically acceptable wormhole solutions that certainly do not violate the weak energy condition ($\rho_0\geq 0$, $\rho_0+p_r\geq 0$, $\rho_0+p_z\geq 0$) and the dominant energy condition ($\rho_0\geq 0$, $-\rho_0\leq p_r\leq \rho_0$, $-\rho_0\leq p_z\leq \rho_0$). To satisfy the requirements of the strong energy condition ($\rho_0\geq 0$, $\rho_0+p_r\geq 0$, $\rho_0+p_z\geq 0$, $\rho_0+2p_r+p_z\geq 0$), the EoS parameters must obey the constraints~\eqref{c1} if $\Lambda<0$ and the following constraints if $\Lambda>0$:
\begin{align}
\label{c5}&\omega_r=-1\;\text{ and }\;\omega_z=-1,\;\text{ or}\\
\label{c6}&-1< \omega_r< -\frac{1}{5}\;\text{ and }\;-2\omega_r-1\leq\omega_z\leq 1,\;\text{ or}\\
\label{c7}&-\frac{1}{5}\leq\omega_r\leq\frac{1}{3}\;\text{ and }\;\frac{\omega_r-1}{2}<\omega_z\leq 1,\;\text{ or}\\
\label{c8}&\frac{1}{3}<\omega_r<1\;\text{ and }\;2\omega_r-1<\omega_z\leq 1.
\end{align}

Considered as interiors all the matter should be distributed inside cylindrical surfaces $\Sigma_+$ and $\Sigma_-$ of finite radii $u_+>u_0$ and $u_-<-u_0$. Referring to~\cite{Bonnor,Whittaker}, the effective mass $m_+$ and angular momentum $j_+$ per unit $z$-coordinate length of matter enclosed by $\Sigma_+$ are defined by (there are similar definitions for $m_-$ and $j_-$ concerning the matter distribution enclosed by $\Sigma_-$):
\begin{align}
&m_+=2\pi \int_{0}^{u_+}\tau_{t}{}^\mu n_\mu\sqrt{|g|}~{\rm d}u,\\
&j_+=-2\pi \int_{0}^{u_+}\tau_{\phi}{}^\mu n_\mu\sqrt{|g|}~{\rm d}u,
\end{align}
where $\tau_{\nu}{}^\mu$ is the total SET given in~\eqref{tau} and $g$ is the determinant of the metric. The vector $n_\mu$ is the unit normal to the spacelike surface of integration, which is the hypersurface $t=\text{const.}$ yielding $n_\mu=\delta_{\mu}^t/\sqrt{g^{tt}}$ with $g^{tt}$ being the component $tt$ of the inverse metric. Using~\eqref{metric2} along with $g=-1/q^2$ we find
\begin{align}
\label{mp}&m_+=\frac{2\pi u_0}{q} \Big(\rho_0-\frac{\Lambda}{\kappa}\Big) \int_{0}^{x_+}\sqrt{\frac{x^2+1}{\delta x^2+1}}~{\rm d}x,\\
\label{jp}&j_+=\frac{4\pi u_0^2\rho_0\omega_0(1+\omega_r)}{q^2} \int_{0}^{x_+}x\sqrt{\frac{x^2+1}{\delta x^2+1}}~{\rm d}x,
\end{align}
where $p_\phi=p_r=\omega_r \rho_0$ has been used, $x\equiv u/u_0$ and $x_+\equiv u_+/u_0$. These integrals {are expressed in terms of the incomplete elliptic integral of the second kind and the log function as given in the Appendix}. It is obvious that $j_+$ has the sign of $\omega_0$ and is zero in the static case (no rotation) and in the extreme case $\omega_r=-1$. The cosmological constant being small, we expect that in physically interesting situations to have $\rho_0>|\Lambda|/\kappa$ and thus $m_+>0$.

The constants of integration ($\omega_0,\,u_0/q$) and $\rho_0\propto \omega_0^2$~\eqref{o02} {could be determined analytically in terms of ($m_+,\,j_+,\,x_+$) if the expressions of ($m_+,\,j_+$), as given in the Appendix, were inversible. This is, however, possible for $x_+\ll 1$ where ($m_+,\,j_+$) expand as (for all $\delta$):
\begin{align*}
&m_+=\frac{\pi  u_0 (\kappa  \rho _0-\Lambda ) x_+ [6+(1-\delta ) x_+^2]}{3 q \kappa }+\mathcal{O}(x_+^5),\\
&j_+=\frac{\pi  u_0^2 \rho _0 \omega _0 (1+\omega _r) x_+^2 [4+(1-\delta ) x_+^2]}{2 q^2}+\mathcal{O}(x_+^6).
\end{align*}
On eliminating the ratio $u_0/q$ and using $\rho_0=2\omega_0^2/[\kappa(1+\omega_{z})]$~\eqref{o02}, we arrive at
\begin{multline*}
9 \kappa  m_+^2  (1+\omega _z) (1+\omega _r)[4+(1-\delta ) x_+^2] \omega _0^3\\ = \pi  j_+ [6+(1-\delta ) x_+^2]^2
[\Lambda  (1+\omega _z)-2 \omega _0^2]^2,
\end{multline*}
which can be solved by radicals for $\omega_0$ in terms of ($m_+,\,j_+,\,x_+$). The roots of this quartic equation in $\omega_0$, which can be obtained using a computer algebra system, are sizable and cannot be reproduced in this work. Then, the ratio $u_0/q$ is obtained from the above-given expansion of $m_+$. Finally, using the expression of $q$~\eqref{form} we can determine all the constants of integration in terms of the physical quantities ($m_+,\,j_+,\,x_+,\,\omega_{r},\,\omega_{z},\,\Lambda,\,\kappa$)}.

The curvature scalar, $\mathcal{R}$, and the Kretschmann scalar, $R^{\mu\nu\alpha\beta}R_{\mu\nu\alpha\beta}$, assume the following expressions:
\begin{equation}\label{scalar}
\mathcal{R}=2 (\omega_0^2-q^2),
\end{equation}
\begin{align}
&R^{\mu\nu\alpha\beta}R_{\mu\nu\alpha\beta}=4(q^4-6\omega_0^2q^2+11\omega_0^4),
\end{align}
which are finite constants.

\section{Cosmological rotating solutions with two axes of symmetry at finite proper distance\label{sectas}}
If the constant $q^2$ in~\eqref{form} is negative we set $Q^2=-q^2$ in~\eqref{form} and the solution yields
\begin{equation}\label{alphan}
	{\rm e}^{2 \alpha }={\rm e}^{2 \beta }=\frac{c^2\text{sech}^2(c r)}{Q^2},\quad \Omega =\frac{2c\omega_0}{Q^2}~\tanh (cr),
\end{equation}
where we have dropped an additive constant $\Omega_0$ in the expression of $\Omega$. On setting $u=c\tanh(cr)/Q$ and $u_0^2=c^2/Q^2$, we bring the metric to the form
\begin{multline}\label{metricQ}
	{\rm d}s^2=  \Big[ {\rm d}t -\frac{2c\omega_0\tanh (cr)}{Q^2}~{\rm d}\phi \Big]^2- \frac{c^2\text{sech}^2(c r)}{Q^2}~{\rm d}r^2 \\- \frac{c^2\text{sech}^2(c r)}{Q^2}~{\rm d}\phi^2
	- {\rm d}z^2,\qquad\qquad\qquad\qquad
\end{multline}
\begin{multline}\label{metric2Q}
	{\rm d}s^2=  \Big( {\rm d}t -\frac{2\omega_0 u}{Q}~{\rm d}\phi \Big)^2- \frac{1}{Q^2(u_0^2-u^2)}~{\rm d}u^2 \\- (u_0^2-u^2){\rm d}\phi^2
	- {\rm d}z^2.\qquad\qquad\qquad\qquad
\end{multline}
Note that this solution shares with the wormhole solution~\eqref{metric}-\eqref{metric2} the same values of the physical quantities ($\rho_0,\;\omega_0^2$), given in~\eqref{rho1} and~\eqref{o02}. The corresponding TEGR solution is also given by~\eqref{metricQ}-\eqref{metric2Q} and shares the same value of the torsion $T=2\omega_0^2$ with the wormhole solution~\eqref{metric}-\eqref{metric2}.

The circular radius ${\rm e}^{\beta}=\sqrt{u_0^2-u^2}$ vanishes at $u=\pm u_0$ signaling the presence of two axes of symmetry. In the vicinity of these two axes the solution~\eqref{metricQ}-\eqref{metric2Q} has CTCs where $g_{\phi\phi}=-(u_0^2-u^2-4\omega_0^2u^2/Q^2)$ becomes positive (the nonrotating solution, $\omega_0=0$, has no CTCs). The integral \[\int_{-u_0}^{u_0}\frac{1}{Q\sqrt{u_0^2-u^2}}~{\rm d}u=\frac{\pi}{Q},\] being convergent, the two axes are at finite proper distance from each other and there is no spatial infinity. The absence of spatial infinity is also known for the nonrotating Melvin solution~\cite{exact,exact2,Melvin}, however, for the latter the proper distance of the two axes of symmetry \[\int_{0}^{\infty}(1+B^2u^2){\rm d}u,\]
diverges ($B$ is the magnetic field). Such nontrivial behaviors of the intrinsic geometry are familiar with static and rotating, cylindrically symmetric and/or axially symmetric, metrics and more examples are provided in~\cite{exact}.

\section{Conclusion\label{secc}}
As we mentioned in the Introduction, the determination of rotating and static solutions around an infinite axis is still attracting much attention. In this work, we presented a first set of two cosmological (energy density constant) rotating solutions in GR and TEGR gravity sourced by anisotropic fluids (isotropic in a plane perpendicular to the axis of symmetry) and extended the existing list of solutions pertaining to GR. These solutions have the property that their angular velocity is proportional to the cosmological constant.

We have shown that the EoS parameters obey a large set of values ensuring the satisfaction of all local energy conditions for the rotating wormholes. Such solutions can straightforwardly be matched to exterior rotating Minkowskian metrics.

The other cosmological rotating solution has two axes of symmetry at finite proper distance where one axis can be regularized upon appropriately fixing the value of the additive constant of integration in the expression of $\Omega$ (which we have dropped) at the expense of rendering the time coordinate periodic.

\section*{Appendix: Analytical expression for ($\pmb{m_+,\,j_+}$)\label{seca}}
\renewcommand{\theequation}{A.\arabic{equation}}
\setcounter{equation}{0}
{The analytical expressions of ($m_+,\,j_+$) depend on the sign of $\delta$, however, their series expansions about $x_+=0$ do not, as this is shown in the new paragraph following~\eqref{jp}. For $\delta>0$ we have
\begin{align*}
&m_+=-\frac{2 {\rm i} \pi  u_0 (\kappa  \rho _0-\Lambda )}{\kappa  q \sqrt{\delta }}~ E\Big({\rm i} \ln (\sqrt{\delta } x_++\sqrt{1+\delta
	x_+^2})\big|\frac{1}{\delta }\Big),\\
&j_+=\frac{2 \pi  u_0^2 \rho _0 \omega _0 (1+\omega _r)}{\delta^{3/2}  q^2}\bigg[\sqrt{1+\delta  x_+^2} \sqrt{\delta  (1+x_+^2)}-\sqrt{\delta }\\
&+(1-\delta )  \ln \bigg(\frac{1+\sqrt{\delta
}}{\sqrt{1+\delta  x_+^2}+\sqrt{\delta  (1+x_+^2)}}\bigg)\bigg].
\end{align*}
For $\delta=0$, we have
\begin{align*}
&m_+=\frac{\pi  u_0 (\kappa  \rho _0-\Lambda ) [x_+ \sqrt{1+x_+^2}+\ln (x_++\sqrt{1+x_+^2})]}{\kappa  q},\\
&j_+=\frac{4 \pi  u_0^2 \rho _0 \omega _0 (1+\omega _r) [(1+x_+^2) \sqrt{1+x_+^2}-1]}{3 q^2}.
\end{align*}
For $\delta<0$ we obtain
\begin{align*}
&m_+=\frac{2 \pi  u_0 (\kappa  \rho _0-\Lambda )}{\kappa  q \sqrt{-\delta }}~ E\Big(-{\rm i} \ln ({\rm i}\sqrt{-\delta } x_++\sqrt{1+\delta
		x_+^2})\big|\frac{1}{\delta }\Big),\\
&j_+=\frac{2 \pi  u_0^2 \rho _0 \omega _0 (1+\omega _r)}{\delta^{3/2}  q^2}\bigg[{\rm i}\sqrt{1+\delta  x_+^2} \sqrt{-\delta  (1+x_+^2)}-{\rm i}\sqrt{\delta }\\
&+(1-\delta )  \ln \bigg(\frac{1+{\rm i}\sqrt{-\delta
}}{\sqrt{1+\delta  x_+^2}+{\rm i}\sqrt{-\delta  (1+x_+^2)}}\bigg)\bigg],
\end{align*}
where $0<x_+<1/\sqrt{-\delta}$ if $\delta<0$}.

{In all these expressions ${\rm i}^2=-1$ and \[E(z |m)\equiv\int_{0}^{z}\sqrt{1-m\sin^2\theta}~{\rm d}\theta ,\] is the incomplete elliptic integral of the second kind. Here $z$ is generally a complex number input. Despite a complex input, the output of the above-given expressions of ($m_+,\,j_+$) is always real for all $\delta$ (for $\delta<0$ the output is real provided $0<x_+<1/\sqrt{-\delta}$).}


\section{Introduction}
Quadrature amplitude modulation (QAM) is a family of modulation schemes, which are widely used in today's communication systems due to their spectral efficiencies. Communication systems employing QAM requires a phase recovery for signal demodulation at the receiver. This phase recovery process is often not perfect due to hardware impairments and/or imperfect channel estimation \cite{Gappmair2017}. Consequently, phase noise can potentially be introduced to the received waveform causing an additional performance loss.

Regarding the impact of phase noise on the bit-error rate (BER) performance, there has been a number of studies for $M$-ary phase-shift keying (PSK) modulation ($M$-PSK) \cite{Jang2013,Song2015,Gappmair2017}, and references therein. Remarkably, for the case of rectangular QAM, an exact closed-form conditional BER expression under AWGN channels has been reported very recently \cite{Jafari2020
where the derivation follows the similar approach presented in \cite{Cho2002} for the case without phase noise. Specifically, exact expressions were derived for the first two bits of symbols on the I channel. General expressions for arbitrary bits of the I and Q channels were then obtained based on the regularities of the particular Gray code sequence (labeling) being considered. It should be noted that there is a number of ways to label symbols using the Gray code in a QAM constellation, resulting in different BERs. As shown in \cite{Agrell2004} that under particular assumption on the channel, binary reflected Gray code (BRGC) sequences are the optimal labeling scheme for $M$-PSK, $M$-ary pulse-amplitude modulation ($M$-PAM), and rectangular QAM. Although BRGC sequences are equivalent in terms of the BER, each gives rise to a particular characterization of bit decision regions and positions, resulting in a different closed-form expression. Note that the authors in \cite{Jafari2020,Cho2002} assumed the original BRGC sequence proposed by Frank Gray in \cite{Gray1953} for their derivations.

In this paper, we aim at providing a rigorous derivation for the conditional BER of rectangular QAM over AWGN channels under the presence of phase noise. The obtained expression serves as a benchmark for evaluating the average BER in a future study.  In our previous work \cite{Nguyen2020}, we gave a sketch of the derivation and numerically calculated the BERs for two specific cases of $8 \times 4$ and $4 \times 4$ QAM. This work attempts at a complete analysis for the case of general rectangular QAM. Specifically, we first discuss in Section II the relationship between BRGC sequences and characterizations of bit decision regions in the case of $M$-PAM. This serves as the motivation for choosing the original BRGC sequence in our analysis. Unlike the previous works where closed-form expressions were obtained through generalization from specific cases \cite{Jafari2020,Cho2002}, our approach is to derive an exact closed-form expression for an $M$-PAM
by analytically characterizing the bit decision regions and positions. The resulted expression can then be naturally applied to compute the conditional BER of a rectangular QAM conditioned on a given phase noise in Section III. Numerical and simulation results are given in Section IV to verify the accuracy of the analysis and to illustrate the impact of phase noise on the conditional BER.

\section{Bit Error Rate of PAM}
\subsection{Signal Model}
Let us consider an $M$-ary PAM, where $\log_2M$ is a positive integer.
The signal waveforms can be expressed as
\begin{align}
s(t) = A_M\cos2\pi f_c t, \hspace{5mm} 0 \leq t < T,
\end{align}
where $A_M$ is the signal amplitude, which is selected from the set $\{\pm d, \pm 3d, \hdots, \pm (M- 1)d\}$. $f_c$ is the carrier frequency, and $T$ is the symbol duration \cite{Cho2002}. The $M$ symbols are drawn from a BRGC sequence  with the distance between any two adjacent symbols being $2d$, which is given by
\begin{align}
d = \sqrt{\frac{3\log_2(M)E_b}{M^2-1}},
\end{align}
where $E_b$ is the bit energy. The received waveform after being distorted by an additive noise is then written as
\begin{align}
r(t) = s(t) + n(t),
\end{align}
where $n(t)$ denotes a zero-mean AWGN with two-sided power spectrum density $N_0/2$. In the following, the time index $t$ is omitted for simplicity.
\begin{figure}[h]
\centering
\includegraphics[width = 8.8cm, height = 2.5cm]{Figures/00_GrayCode_a.eps}
\caption{Standard construction of $\mathcal{G}_3$.}
\label{GrayConstruction}
\end{figure}
\begin{figure*}[bt]
\centering
\begin{subfigure}[b]{0.48\textwidth}
\includegraphics[width = 8.2cm, height = 2.8cm]{Figures/01_BitDemaping_a.eps}
\caption{Bit decision regions and positions according to $\mathcal{G}_3$.}
\label{bit_demaping_1}
\end{subfigure}
\centering
\begin{subfigure}[b]{0.48\textwidth}
\includegraphics[width = 8.2cm, height = 2.8cm]{Figures/01_BitDemaping_b.eps}
\caption{Bit decision regions and positions according to $\mathcal{G}'_3$.}
\label{bit_demaping_2}
\end{subfigure}
\caption{An example of two different 8-PAM signal point arrangements.}
\label{decRegion-8PAM}
\end{figure*}
\subsection{Bit Decision Regions}
\subsubsection{Construction of Gray code sequences}
Since our analytical framework directly depends on Gray code sequences' construction, specifically the BRGC, we discuss this issue in this section. For the sake of readability, let us revisit the construction of BRGC sequences proposed in \cite{Gray1953}. According to this approach, let $\mathcal{G}_1 = (1, 0)$ be the trivial 1-bit sequence. For $n \geq 2$, denote $\mathcal{G}_{n-1}$ as an $(n-1)$-bit BRGC sequence and $\overline{\mathcal{G}}_{n-1}$ as the reflected version of $\mathcal{G}_{n-1}$ (i.e., codewords of $\mathcal{G}_{n-1}$ are placed in the reverse order). An extra bit `1' and `0' are added to each codeword of $\mathcal{G}_{n-1}$ and $\overline{\mathcal{G}}_{n-1}$ from the left (i.e., prefixing), respectively. An $n$-bit BRGC sequence $\mathcal{G}_n$ is then obtained by concatenating the two newly formed sequences. We regard this as the standard construction, and the resulted sequence is called the standard BRGC sequence. For example, Fig. \ref{GrayConstruction} illustrates the generation of the standard 3-bit BRGC sequence $\mathcal{G}_3$. Note that, from the standard $\mathcal{G}_n$ sequence, it is possible to generate another $n!2^n - 1$ sequences from it by means of permutation and/or complementation of all the codewords' bits in the same manner \cite{Arazi1984}. For example, permuting the first and third bits and complementing the second bit in all codewords of $\mathcal{G}_3$ result in a new Gray code sequence $\mathcal{G}'_3 = (101~001~011~111~110~010~000~100)$.

In this paper, we are interested in an $M$-PAM where the arrangement of symbols\footnote{We use the terms ``symbols" and ``codewords" interchangeably.} is drawn from the standard $\mathcal{G}_{\log_2(M)}$ BRGC sequence. By convention, codewords in the first and second half of the sequence are placed in the left (negative $A_M$ values) and right (positive $A_M$ values) half plane of the signal space, respectively. Denote $K = \log_2(M)$ and let $b_{K,k}$ be the symbol's $k-$th bit from the left ($k \leq K$). It is worth mentioning that different code sequences lead to different decision regions of $b_{K, k}$ as illustrated in Fig. \ref{decRegion-8PAM} for the case of 8-PAM with symbol arrangements according to $\mathcal{G}_3$ and $\mathcal{G}'_3$. In this example, although the two sequences result in different decision regions of each bit (hence its error rate), the overall BERs of the two cases are the same.  In fact, \cite{Agrell2004} shows that trivial operations on codewords of a sequence such as bit permutation and complementation do not affect its BER. Hence, it suffices to derive the BER of an $M$-PAM whose symbol arrangement is drawn from the standard $\mathcal{G}_{K}$ sequence.
\subsubsection{Analytical characterization of bit decision regions}
To facilitate the discussion, the following definition and proposition are necessary.
\begin{Definition}
Let $b^i_{K, 1}b^i_{K,2}...b^i_{K,K}$ ($b^i_{K, k} \in \{0, 1\}$) be the $i-$th codeword of $\mathcal{G}_K$. Then, $\mathcal{B}_{K, k} = \left(b^1_{K, k}, b^2_{K, k}, \cdots, b^{2^K}_{K, k}\right)$ is defined as a binary $2^K$-tuple  whose $i-$th element  is the value of  $b_{K, k}$ in the $i-$th codeword of $\mathcal{G}_K$.
\end{Definition}
\begin{Proposition}
 \begin{align}
\mathcal{B}_{K, k} \!=\! \left\{ {\begin{array}{*{20}{l}}
\!\!\!(\underbrace{1, \cdots, 1}_{2^{K-1}}, \underbrace{0, \cdots, 0}_{2^{K-1}}) \hspace*{0pt}\hfill {{\text{if~~}}} k = 1, \\
\!\!\!(\underbrace{1, \cdots, 1}_{2^{K-k}}, \underbrace{0, \cdots, 0}_{2^{K-k + 1}}, \underbrace{1, \cdots, 1}_{2^{K-k + 1}}, \cdots, \underbrace{0, \cdots, 0}_{2^{K-k + 1}},  \underbrace{1, \cdots, 1}_{2^{K-k }}) \nonumber \\ \hspace*{0pt}\hfill {{\text{if~~}}}  k \geq 2. \\
\end{array}} \right.
\end{align}
\end{Proposition}

\begin{proof} Due to the construction rule of $\mathcal{G}_K$ from $\mathcal{G}_{K-1}$ which has $2^{K-1}$ codewords, it is obvious to see that $\mathcal{B}_{K, 1} = (\underbrace{1, \cdots, 1}_{2^{K-1}}, \underbrace{0, \cdots, 0}_{2^{K-1}})$. For $k \geq 2$, let $\overline{\mathcal{B}}_{K, k}$ be the reflected version of $\mathcal{B}_{K, k}$. Notice that $\mathcal{B}_{K, k}$ is symmetrical since $\mathcal{B}_{K, k}  = \left(\mathcal{B}_{K-1, k-1}\overline{\mathcal{B}}_{K-1, k-1}\right)$. Thus, $\mathcal{B}_{K, k} = \overline{\mathcal{B}}_{K, k}$ and we can write
\begin{align}
\mathcal{B}_{K, k} & = \left(\mathcal{B}_{K-1, k-1}\overline{\mathcal{B}}_{K-1, k-1}\right) \nonumber \\
			     & = \left(\mathcal{B}_{K-2, k-2}\overline{\mathcal{B}}_{K-2, k-2}{\mathcal{B}}_{K-2, k-2}\overline{\mathcal{B}}_{K-2, k-2}\right) \nonumber \\
			     & ~~\vdots \nonumber \\
			     & = (\underbrace{\mathcal{B}_{K-k+1, 1}\overline{\mathcal{B}}_{K-k+1, 1}\cdots\mathcal{B}_{K-k+1,1}\overline{\mathcal{B}}_{K-k+1, 1}}_{2^{k-1} \text{pairs}~\mathcal{B}_{K-k+1,1}\overline{\mathcal{B}}_{K-k+1, 1}} ).
\label{bit-pattern}
\end{align}
Since $\mathcal{B}_{K-k+1, 1} = (\underbrace{1, \cdots, 1}_{2^{K-k}}, \underbrace{0, \cdots, 0}_{2^{K-k}})$ and $\overline{\mathcal{B}}_{K-k+1, 1} = (\underbrace{0, \cdots, 0}_{2^{K-k}}, \underbrace{1, \cdots, 1}_{2^{K-k}})$, the preposition is proved following the observation from \eqref{bit-pattern}. \qedhere
\end{proof}

The proposition shows in the trivial case of $k = 1$ that $b_{K, 1} = 1$ if $r \leq 0$ and $b_{K, 1} = 0$ if $r > 0$. For $k \geq 2$, a decision rule for $b_{K, k}$ can be written as
\begin{align}
b_{K, k} = \left\{ {\begin{array}{*{20}{c}}
1 & {{\text{if~~}}}& r \odot \mathcal{R}^{(1)}_{K, k}, \\
0 & {{\text{if~~}}}& r \odot \mathcal{R}^{(0)}_{K, k}, \\
\end{array}} \right.
\end{align}
where $\mathcal{R}^{(1)}_{K, k} = \left\{\left(\delta^{(1), l}_{K, k}d,~\delta^{(1), u}_{K, k}d\right)~\Big|~\delta^{(1), l}_{K, k} \in \left\{\mathbb{Z}, -\infty\right\}, \right.$ $\left. \delta^{(1), u}_{K, k} \in \left\{\mathbb{Z}, +\infty\right\}\right\}$ and $\mathcal{R}^{(0)}_{K, k} = \left\{\left(\delta^{(0), l}_{K, k}d,~\delta^{(0), u}_{K, k}d\right)~\Big|~\delta^{(0), l}_{K k} \in \mathbb{Z},  \delta^{(0), u}_{K, k} \in \mathbb{Z}\right\}$ are the sets of decision regions of bits `1' and `0', respectively. Note that, bit `1's decision regions at the two extremes of the signal space are unbounded. This is represented by $ \delta^{(1), l}_{K, k} = -\infty$ and $ \delta^{(1), u}_{K, k} = +\infty$. With $k = 3$ from the above example, $\mathcal{R}^{(1)}_{3, 3} = \left\{\left(-\infty, -6d\right), ~\left(-2d, 2d\right),~\left(6d, +\infty\right)\right\}$ and $\mathcal{R}^{(0)}_{3, 3} = \left\{\left(-6d, -2d\right), ~\left(2d, 6d\right)\right\}$. We use the notation $r  \odot \mathcal{R}^{(1)}_{K, k}$ (or  $r  \odot \mathcal{R}^{(0)}_{K, k}$) to imply that there exist some $\left(\delta^{(1), l}_{K, k}d,~\delta^{(1), u}_{K, k}d\right) \in \mathcal{R}^{(1)}_{K, k}$ that satisfies $\delta^{(1), l}_{K, k}d \leq r <  \delta^{(1), u}_{K, k}d$. We now derive analytical expressions of $\left\{\delta^{(0), l}_{K, k}\right\}$,  $\left\{\delta^{(0), u}_{K, k}\right\}$, $\left\{\delta^{(1), l}_{K, k}\right\}$, and $\left\{\delta^{(1), u}_{K, k}\right\}$.

A direct implication of the proposition is that the number of `1's and `0's regions of $b_{K, k}$ is $\left|\mathcal{R}^{(1)}_{K, k}\right| + \left|\mathcal{R}^{(0)}_{K, k}\right| = 2^{k-1}+ 1$.
Accordingly, let $\mathbf{p}_k = \begin{bmatrix}-1 & 0 & 1 & \cdots & 2^{k-1}-1 & 2^{k-1} \end{bmatrix}$ be the vector whose two consecutive elements represent the starting and ending points of a region.
Denote $\Delta(x, K, k)$ as a function of $x$, $K$, and $k$, which returns the values of $\left\{\delta^{(0), l}_{K, k}\right\}$,  $\left\{\delta^{(0), u}_{K, k}\right\}$, $\left\{\delta^{(1), l}_{K, k}\right\}$, and $\left\{\delta^{(1), u}_{K, k}\right\}$ for $x$ taking on values from $\mathbf{p}_k$.
To account for $ \delta^{(1), l}_{K, k} = -\infty$ and $ \delta^{(1), u}_{K, k} = +\infty$, it is  mathematically convenient to use the convention $\Delta(-1, K, k) = -\infty$ and $\Delta(2^{k-1}, K, k) = +\infty$.
Observe that there are $2^{K-k}$ bits `1' in the first decision regions followed by interleaving regions of `0' and `1' each containing $2^{K-k + 1}$ bits  and ended by a region (`0' if $k = 1$ or `1' if $k \geq 2$) having $2^{K-k}$ bit. As the distance between any two adjacent codewords is $2d$ with the first one starting at $(-2^K+1)d$, the observation gives rise to the definition of $\Delta(x, K, k)$ as
\begin{align}
&\Delta(x, K, k) = \nonumber \\
& \left\{ {\begin{array}{*{5}{l}}
\!\!\!-\infty & , x = -1, \\
\!\!\!-2^K\left(1 - \frac{1}{2^{k-1}}\right) + x\times2^{K - k +2} &, -1 < x < 2^{k-1} , \\
\!\!\!+\infty &, x =  2^{k-1}.
\end{array}} \right.
\label{delta}
\end{align}
Then, $\left\{\delta^{(1), l}_{K, k}\right\}$, $\left\{\delta^{(1), u}_{K, k}\right\}$, $\left\{\delta^{(0), l}_{K, k}\right\}$, and  $\left\{\delta^{(0), u}_{K, k}\right\}$ are given by
\begin{align}
&\left\{\delta^{(1), l}_{K, k}\right\} = \Delta\left(2\mathbf{p}_k^{(1)}+1, K, k\right),
\label{delta_1_l}\\
&\left\{\delta^{(1), u}_{K, k}\right\} = \Delta\left(2\mathbf{p}_k^{(1)}+2, K, k\right),
\label{delta_1_u}\\
&\left\{\delta^{(0), l}_{K, k}\right\} = \Delta\left(2\mathbf{p}_k^{(0)} + 2, K, k\right),
\label{delta_0_l}\\
&\left\{\delta^{(0), u}_{K, k}\right\} = \Delta\left(2\mathbf{p}_k^{(0)}+3, K, k\right),
\label{delta_0_u}
\end{align}
where
\begin{align}
\mathbf{p}_k^{(1)} = \left\{ {\begin{array}{*{20}{l}}
-1 & {{\text{if~~}}}  & k = 1,\\
\begin{bmatrix}-1 & 0 & \cdots  & 2^{k-2} - 1\end{bmatrix} & {{\text{if~~}}} & k \geq 2,
\end{array}} \right.
\label{region-1-index}
\end{align}
and
\begin{align}
\mathbf{p}_k^{(0)} = \left\{ {\begin{array}{*{20}{l}}
-1 & {{\text{if~~}}}  & k = 1,\\
\begin{bmatrix}-1 & 0 & \cdots & 2^{k-2} - 2 \end{bmatrix}& {{\text{if~~}}} & k \geq 2.
\end{array}} \right.
\label{region-0-index}
\end{align}
\subsection{Bit Positions and Conditional BER Expression}
\begin{figure*}[ht]
\begin{align}
&P_{\text{PAM}}(b_{K, k} )   = \text{Pr}\left(r \odot \mathcal{R}^{(1)}_{K, k} \Big| s \in \mathcal{A}^{(0)}_{K, k} \right)
               +  \text{Pr}\left(r \odot \mathcal{R}^{(0)}_{K, k} \Big| s \in \mathcal{A}^{(1)}_{K, k} \right)
               = \frac{1}{{MK}}\left(\sum_{s \in \mathcal{A}^{(0)}_{K, k} }\!\text{Pr}\left(r \odot \mathcal{R}^{(1)}_{K, k} \right)
               \! + \!\sum_{s \in \mathcal{A}^{(1)}_{K, k} }\!\text{Pr}\left(r \odot \mathcal{R}^{(0)}_{K, k} \right) \right) \nonumber \\
              & = \frac{1}{MK}\left(\sum_{\overline{p}^{(0)} \in \mathbf{\overline{p}}_k^{(0)}}\sum_{j^{(0)} \in \mathbf{j}^{\overline{p}^{(0)}}_{K, k}}\sum_{p^{(1)} \in \mathbf{p}_k^{(1)}}\text{Pr}\left(\delta^{p^{(1)}, l}_{K, k}d \leq A^{\overline{p}^{(0)}, j^{(0)}}_{K, k}d + n < \delta^{p^{(1)}, u}_{K, k}d\right) \right. \nonumber \\
              &~~~~~~~~~ \left.+ \sum_{\overline{p}^{(1)} \in \mathbf{\overline{p}}_k^{(1)}}\sum_{j^{(1)} \in \mathbf{j}^{\overline{p}^{(1)}}_{K, k}}\sum_{p^{(0)} \in \mathbf{p}_k^{(0)}}\text{Pr}\left(\delta^{p^{(0)}, l}_{K, k}d \leq A^{\overline{p}^{(1)}, j^{(1)}}_{K, k}d + n < \delta^{p^{(0)}, u}_{K, k}d\right) \right)  \nonumber \\
              & = \frac{1}{2MK}\left(\sum_{\overline{p}^{(0)} \in \mathbf{\overline{p}}_k^{(0)}}\sum_{j^{(0)} \in \mathbf{j}^{\overline{p}^{(0)}}_{K, k}}\sum_{p^{(1)} \in \mathbf{p}^{(1)}_k}\text{erfc}\left(\frac{d}{\sqrt{N_0}}\left(\delta^{p^{(1)}, l}_{K, k} - A^{\overline{p}^{(0)}, j^{(0)}}_{K, k}\right)\right) - \text{erfc}\left(\frac{d}{\sqrt{N_0}}\left(\delta^{p^{(1)},u}_{K, k} - A^{\overline{p}^{(0)}, j^{(0)}}_{K, k}\right)\right) \right . \nonumber \\ & \left. ~~~~~~~~~~~ +\sum_{\overline{p}^{(1)} \in \mathbf{\overline{p}}_k^{(1)}}\sum_{j^{(1)} \in \mathbf{j}^{\overline{p}^{(1)}}_{K, k}}\sum_{p^{(0)} \in \mathbf{p}^{(0)}_k}  \text{erfc}\left(\frac{d}{\sqrt{N_0}}\left(\delta^{p^{(0)}, l}_{K, k} - A^{\overline{p}^{(1)}, j^{(1)}}_{K, k} \right)\right)  - \text{erfc}\left(\frac{d}{\sqrt{N_0}}\left(\delta^{p^{(0)}, u}_{K, k} - A^{\overline{p}^{(1)}, j^{(1)}}_{K, k}\right)\right)  \vphantom{\sum_{\left(\delta^{(1), l}_{K, k}d, \delta^{(1), u}_{K, k}d\right) \in \mathcal{R}^{(1)}_{K, k} }\sum_{q = -\frac{M}{2} }^{\frac{M}{2} - 1}} \right).
\label{BER-PAM}
\end{align}
\noindent\makebox[\linewidth]{\rule{\textwidth}{1pt}}
\end{figure*}
To formulate an expression for the BER of $b_{K, k}$, it is now necessary to determine the positions of the elements of $\mathcal{B}_{K, k}$ in the signal space.
In order to so, let $\mathcal{A}^{(1)}_{K, k} = \left\{A^{(1)}_{K, k}d ~ \big| ~ A^{(1)}_{K, k} \in \mathbb{Z} \right\}$ and $\mathcal{A}^{(0)}_{K, k} = \left\{A^{(0)}_{K, k}d ~ \big| ~ A^{(0)}_{K, k} \in \mathbb{Z} \right\}$ be the sets of positions of the symbols whose $b_{K, k}$'s are `1' and `0', respectively. In our example of $8$-PAM in Fig. \ref{bit_demaping_1}, for $k = 3$, we have $\mathcal{A}^{(1)}_{3, 3} = \left\{-7d, -d, d, 7d  \right\}$ and $\mathcal{A}^{(0)}_{3, 3} = \left\{-5d, -3d, 3d, 5d  \right\}$. Since the values of $\mathcal{A}^{(1)}_{K, k}$ $(\mathcal{A}^{(0)}_{K, k})$ lie on $\mathcal{R}^{(1)}_{K, k}$ $(\mathcal{R}^{(0)}_{K, k})$, they can be determined using the defined $\Delta(\cdot)$ function. To make characterizations of  $\mathcal{A}^{(0)}_{K, k}$ and $\mathcal{A}^{(1)}_{K, k}$ convenient,
we bound the two unbounded regions at $-2^K$ and $2^K$, which correspond to $\Delta\left(-\frac{1}{2}, K, k\right)$ and $\Delta\left(2^{k-1} - \frac{1}{2}, K, k\right)$, respectively. This naturally necessitates a definition of $\mathbf{\overline{p}}_k = \begin{bmatrix}-\frac{1}{2} & 0 & 1 &\cdots & 2^{k-1} - 1& 2^{k-1} - \frac{1}{2} \end{bmatrix}$.
Hence, we can write $\mathcal{A}^{(1)}_{K, k} = \left\{\left\{A_{K, k}^{(1)}\right\}_{\left[\mathbf{\overline{p}}^{(1)}_k\right]_i}\right\}$ (similarly for $\mathcal{A}^{(0)}_{K, k}$), where $\left\{A_{K, k}^{(1)}\right\}_{\left[\mathbf{\overline{p}}^{(1)}_k\right]_i}$ is the set of bit `1's positions corresponding to $\left[\mathbf{\overline{p}}^{(1)}_k\right]_i$, which is the $i-$th element of $\mathbf{\overline{p}}^{(1)}_k$ and is defined as
\begin{align}
\mathbf{\overline{p}}_k^{(1)} = \left\{ {\begin{array}{*{20}{l}}
-\frac{3}{4} & {{\text{if~~}}}  & k = 1,\\
\begin{bmatrix}-\frac{3}{4} & 0 & \cdots  & 2^{k-2} - 1\end{bmatrix} & {{\text{if~~}}} & k \geq 2.
\end{array}} \right.
\label{region-1-index}
\end{align}
Similarly, denote $\left[\mathbf{\overline{p}}^{(0)}_k\right]_i$ as the $i-$th element of $\mathbf{\overline{p}}^{(0)}_k$, which is defined to be the same as  $\mathbf{{p}}^{(0)}_k$ in \eqref{region-0-index}.
We also define the length of a region $\left(\Delta\left([\mathbf{\overline{p}}_k]_i, K, k\right)~\Delta\left([\mathbf{\overline{p}}_k]_{i+1}, K, k\right)\right)$ straightforwardly as
\begin{align}
l\left([\mathbf{\overline{p}}_k]_i, K, k\right) &= \Delta([\mathbf{\overline{p}}_k]_{i + 1}, K, k) - \Delta\left([\mathbf{\overline{p}}_k]_i, K, k\right)\nonumber \\
& = \left([\mathbf{\overline{p}}_k]_{i+1} - [\mathbf{\overline{p}}_k]_i\right)\times2^{K - k +2}
\label{delta-modified}
\end{align}
Now, $\left\{A_{K, k}^{(1)}\right\}_{\left[\mathbf{\overline{p}}^{(1)}_k\right]_i} $ and $\left\{A_{K, k}^{(0)}\right\}_{\left[\mathbf{\overline{p}}^{(0)}_k\right]_i}$ are respectively given by
\begin{align}
\left\{A_{K, k}^{(1)}\right\}_{\left[\mathbf{\overline{p}}^{(1)}_k\right]_i}  = \Delta\left(2\left[\mathbf{\overline{p}}^{(1)}_k\right]_i + 1, K, k\right) + 2\mathbf{j}_{K, k}^{\left[\mathbf{\overline{p}}^{(1)}_k\right]_i} + 1,
\label{position_1}
\end{align}
and
\begin{align}
\left\{A_{K, k}^{(0)}\right\}_{\left[\mathbf{\overline{p}}^{(0)}_k\right]_i} = \Delta\left(2\left[\mathbf{\overline{p}}^{(0)}_k\right]_i + 2, K, k\right) + 2\mathbf{j}_{K, k}^{\left[\mathbf{\overline{p}}^{(0)}_k\right]_i} + 1,
\label{position_0}
\end{align}
, where $\mathbf{j}^{\left[\mathbf{\overline{p}}^{(1)}_k\right]_i}_{K, k} =  \begin{bmatrix} 0 & 1 & \cdots & \frac{l\left(2\left[\mathbf{\overline{p}}^{(1)}_k\right]_i + 1, K, k\right)}{2} -1 \end{bmatrix}$ and $\mathbf{j}^{\left[\mathbf{\overline{p}}^{(0)}_k\right]_i}_{K, k} = \begin{bmatrix} 0 & 1 & \cdots & \frac{l\left(2\left[\mathbf{\overline{p}}^{(0)}_k\right]_i + 2, K, k\right)}{2} -1 \end{bmatrix}$.

For each transmitted symbol $s$, $b_{K,k}$ is decoded wrongly when $r$ falls into bit `1's decision regions given that $s$ is chosen from the sets of `0's positions or vice-versa. For expressional convenience, let $\delta^{p^{(1)}, l}_{K, k}$,  $\delta^{p^{(1)}, u}_{K, k}$, $\delta^{p^{(0)}, l}_{K, k}$, and $\delta^{p^{(0)}, u}_{K, k}$ be the elements of $\left\{\delta^{(1), l}_{K, k}\right\}$, $\left\{\delta^{(1), u}_{K, k}\right\}$, $\left\{\delta^{(0), l}_{K, k}\right\}$,  and $\left\{\delta^{(0), u}_{K, k}\right\}$ in \eqref{delta_1_l}-\eqref{delta_0_u}, which correspond to some $p^{(1)}_k \in \mathbf{p}^{(1)}_k$ and $p^{(0)}_k \in \mathbf{p}^{(0)}_k$.
Also, let $A^{\overline{p}^{(1)}, j^{(1)}}_{K, k}$ and $A^{\overline{p}^{(0)}, j^{(0)}}_{K, k}$ are the elements of $\left\{A_{K, k}^{(1)}\right\}_{\left[\mathbf{\overline{p}}^{(1)}_k\right]_i} $ and $\left\{A_{K, k}^{(0)}\right\}_{\left[\mathbf{\overline{p}}^{(0)}_k\right]_i}$, which correspond to some $\overline{p}^{(1)}_{K, k} \in \overline{\mathbf{p}}^{(1)}_{K, k}$, $\overline{p}^{(0)}_{K, k} \in \overline{\mathbf{p}}^{(0)}_{K, k}$, $j^{(1)} \in \mathbf{j}^{\overline{p}^{(1)}_{K, k}}_{K, k}$, and  $j^{(0)} \in \mathbf{j}^{\overline{p}^{(0)}_{K, k}}_{K, k}$. The error probability of $b_{K, k}$ can then be expressed in \eqref{BER-PAM}, which is on top of the previous page. Note that, $\text{erfc}(x) = \frac{2}{\sqrt{\pi}}\int_{x}^{\infty}e^{-t^2}dt$ is the complementary error function.

\section{Bit Error Rate of QAM with Phase Noise}
We now consider a rectangular $M_I \times M_Q$ QAM, whose signal constellation can be seen as a combination of two PAM signalings. The signal waveform is expressed by
\begin{align}
s(t) = A_{I}\cos\left(2\pi f_c t\right) - A_{Q}\sin\left(2\pi f_c t\right),
\end{align}
where  $A_{I}$ and $A_{Q}$ are the signal amplitudes with respect to the I and Q channels, respectively. Also, let $d_I$ and $d_Q$ be the half distances between two adjacent symbols on the two channels. Without lost of generality, we assume that $d_I = d_Q = d$, which is given by
\begin{align}
d = \sqrt{\frac{3\log_2\left(M_IM_Q\right)E_b}{M_I^2 + M_Q^2 - 2}}.
\end{align}
The received waveform after being distorted by an AWGN and phase noise can be written as
\begin{align}
r(t) \!=\! A_{I}\cos\left(2\pi f_c t + \theta_e \right) \!- \! A_{Q}\sin\left(2\pi f_c t + \theta_e \right) + n,
\end{align}
where the phase noise component $\theta_e \in [-\pi ~ \pi]$ caused by the imperfections of the phase-locked loop and/or channel estimation can be modeled by a Tikhonov distribution \cite{Gappmair2017}. The I and Q components of the received waveform are then given by
\begin{align}
r_I&= A_I\cos(\theta_e) - A_Q\sin(\theta_e) + n,   \\
r_Q &= A_I\sin(\theta_e) + A_Q\cos(\theta_e) + n.
\end{align}
Denote $m_I = \log_2(M_I)$ and $m_Q = \log_2(M_Q)$, the conditional BER of the $M_I \times M_Q$ QAM can be calculated  by summing up the conditional BERs of the corresponding $m_I$- and $m_Q$-PAM, resulting in
\begin{align}
P_{\text{QAM}}(\theta_e) = \sum_{k = 1}^{m_I} P_{\text{e}}(i_{m_I, k}|\theta_e) + \sum_{k = 1}^{m_Q} P_{\text{e}}(q_{m_Q,k}|\theta_e),
\end{align}
where $i_{m_I, k}$ and $q_{m_Q, k}$ are the $k-$th bits of the symbols of the $M_I$- and $M_Q$-PAM.
Following the BER expression derived for the $M$-PAM in \eqref{BER-PAM}, the conditional BERs of  $i_{m_I, k}$ and $q_{m_Q, k}$ conditioned on $\theta_e$ denoted as $P_{\text{e}}(i_{m_I, k}|\theta_e)$ and $P_{\text{e}}(q_{m_Q, k}|\theta_e)$ are given by
\begin{align}
&P_{\text{e}}(i_{m_I, k} | \theta_e) = \frac{1}{2(m_I + m_Q)M_I M_Q} \nonumber \\
& \times \left(\sum_{\overline{p}^{(0)} \in \mathbf{\overline{p}}_k^{(0)}}\sum_{j^{(0)} \in \mathbf{j}^{\overline{p}^{(0)}}_{m_I, k}}\sum_{p^{(1)} \in \mathbf{p}^{(1)}_k}\sum_{q = -\frac{M_Q}{2}}^{\frac{M_Q}{2}-1}\right. \nonumber \\ & \left.  ~~~~~~~\Psi\left(\frac{d}{\sqrt{N_0}}, \delta^{p^{(1)}, l}_{m_I, k}, -A^{\overline{p}^{(0)}, j^{(0)}}_{m_I, k}, 2q+1, \theta_e\right) \right. \nonumber \\ & \left. ~~~~~- \Psi\left(\frac{d}{\sqrt{N_0}},\delta^{p^{(1)},u}_{m_I, k}, -A^{\overline{p}^{(0)}, j^{(0)}}_{m_I, k}, 2q+1, \theta_e\right)  \right. \nonumber \\ & \left. ~~+\sum_{\overline{p}^{(1)} \in \mathbf{\overline{p}}_k^{(1)}}\sum_{j^{(1)} \in \mathbf{j}^{\overline{p}^{(1)}}_{m_I, k}}\sum_{p^{(0)} \in \mathbf{p}^{(0)}_k}\sum_{q = -\frac{M_Q}{2}}^{\frac{M_Q}{2}-1} \right. \nonumber \\ & \left.  ~~~~~~~\Psi\left(\frac{d}{\sqrt{N_0}}, \delta^{p^{(0)}, l}_{m_I, k},  -A^{\overline{p}^{(1)}, j^{(1)}}_{m_I, k}, 2q+1, \theta_e \right)  \right. \nonumber \\ & \left. ~~~~~- \Psi\left(\frac{d}{\sqrt{N_0}}, \delta^{p^{(0)}, u}_{m_I, k}, -A^{\overline{p}^{(1)}, j^{(1)}}_{m_I, k}, 2q+1, \theta_e\right)  \vphantom{\sum_{\left(\delta^{(1), l}_{K, k}d, \delta^{(1), u}_{m_I, k}d\right) \in \mathcal{R}^{(1)}_{K, k} }\sum_{q = -\frac{M}{2} }^{\frac{M}{2} - 1}} \!\!\! \right),
\label{BER-QAM-I}
\end{align}
and
\begin{align}
&P_{\text{e}}(q_{m_Q, k} |\theta_e)  = \frac{1}{2(m_I + m_Q)M_IM_Q} \nonumber \\ & \times \left(\sum_{\overline{p}^{(0)} \in \mathbf{\overline{p}}_k^{(0)}}\sum_{j^{(0)} \in \mathbf{j}^{\overline{p}^{(0)}}_{m_Q, k}}\sum_{p^{(1)} \in \mathbf{p}^{(1)}_k}\sum_{i = -\frac{M_I}{2}}^{\frac{M_I}{2}-1} \right. \nonumber \\ & \left. ~~~~~~~\Psi\left(\frac{d}{\sqrt{N_0}}, \delta^{p^{(1)}, l}_{m_Q, k}, -A^{\overline{p}^{(0)}, j^{(0)}}_{m_Q, k}, - (2i+1), \theta_e\right)  \right. \nonumber \\ & \left. ~~~~~- \Psi\left(\frac{d}{\sqrt{N_0}}, \delta^{p^{(1)},u}_{m_Q, k}, -A^{\overline{p}^{(0)}, j^{(0)}}_{m_Q, k}, - (2i +1), \theta_e\right) \right . \nonumber \\ & \left.~~  +\sum_{\overline{p}^{(1)} \in \mathbf{\overline{p}}_k^{(1)}}\sum_{j^{(1)} \in \mathbf{j}^{\overline{p}^{(1)}}_{m_Q, k}}\sum_{p^{(0)} \in \mathbf{p}^{(0)}_k}\sum_{i = -\frac{M_I}{2}}^{\frac{M_I}{2}-1} \right. \nonumber \\ & \left. ~~~~~~~\Psi\left(\frac{d}{\sqrt{N_0}}, \delta^{p^{(0)}, l}_{m_Q, k}, -A^{\overline{p}^{(1)}, j^{(1)}}_{m_Q, k}, - (2i+1), \theta_e\right) \right. \nonumber \\ & \left.  ~~~~~- \Psi\left(\frac{d}{\sqrt{N_0}}, \delta^{p^{(0)}, u}_{m_Q, k}, -A^{\overline{p}^{(1)}, j^{(1)}}_{m_Q, k}, -(2i+1), \theta_e\right)  \vphantom{\sum_{\left(\delta^{(1), l}_{K, k}d, \delta^{(1), u}_{m_Q, k}d\right) \in \mathcal{R}^{(1)}_{K, k} }\sum_{i = -\frac{M}{2} }^{\frac{M}{2} - 1}} \!\!\! \right),
\label{BER-QAM-Q}
\end{align}
respectively. In the above expressions, we define $\Psi(\alpha, x, y, z, \varphi) = \text{erfc}\left(\alpha\left(x+y\cos(\varphi)+z\sin(\varphi)\right)\right)$.

\section{Numerical Examples and Dicussions}
Firstly, we show in Fig. \ref{square_QAM} the conditional BER of square $M$-QAM with different constellation sizes as a function of ${E_b}/{N_0}$. In case of the performance with phase noise, $\theta_e = \pi/180$ is chosen. Similar to the case of PSK, we observed significant performance losses due to phase noise when the constellation size is large. At conditional BER = $10^{-3}$, performance losses are negligible in the case of $4$- and $16$-QAM. However, these are about $0.25$ dB, $1$ dB, and $3.6$ dB in the case of $64$-, $256$-, and $1024$-QAM, respectively.
\begin{figure}[ht]
\centering
\includegraphics[scale = 0.315]{Figures/06_BER_vs_Eb_N0_different_constel_size.eps}
\caption{BER of square QAM for different constellation sizes. }
\label{square_QAM}
\end{figure}

It is well-known that square QAM outperforms other rectangular constellations of the same size (i.e., the product $M_IM_Q$ is kept unchanged) in terms of the BER. Motivated by this, in Fig. \ref{rectangular_QAM}, we examine the impact of phase noise on different rectangular shapes of the same size. In this example, $32\times2$, $16\times4$, and $8\times8$ QAM are taken for demonstration. The result shows that the more balance between signal point sizes on the I and Q channels is, the less performance degradation due to phase noise becomes. For instance, at conditional BER = $10^{-4}$, the performance losses of $32\times2$, $16\times4$, and $8\times8$ QAM are $3.4$ dB, $0.75$ dB, and $0.3$ dB, respectively.
As such, it is evident that the square QAM also suffers less phase noise-induced performance loss than other rectangular configurations.
\begin{figure}[ht]
\centering
\includegraphics[scale = 0.315]{Figures/07_BER_vs_Eb_N0_compare_rec_sq.eps}
\caption{BER of different rectangular QAM  configurations. }
\label{rectangular_QAM}
\end{figure}
\section{Conclusion}
We have derived an exact closed-form expression for the conditional BER of rectangular QAM given a value of phase noise. The derivation was based on analytical characterizations of bit decision regions and positions. Numerical results showed that significant performance losses are possible in the case of large constellation sizes. Besides, with the same constellation size, the square QAM is preferable in terms of both performance and loss due to phase noise.
\bibliographystyle{ieeetr}