\section{Introduction}
\label{sect:intro}

Directional wireless communications over millimeter-wave (mmWave, $30$\,GHz--$300$\,GHz) bands is one of the most significant novelties introduced in the fifth-generation (5G) wireless networks. Adopted by 3rd Generation Partnership Project (3GPP) for New Radio (NR) technology, mmWave communications at the frequencies between $24.25$\,GHz and $52.6$\,GHz allow for the data exchange at the rates of \emph{several gigabits per second (Gbps)}. The initial phase of the standardization process for mmWave NR was finished in 2018, and the first deployments already started to appear in Europe, Japan, the United Kingdom, and the United States.

Along this timeline, the utilization of mmWave communications at $60$\,GHz has also been ratified by the Institute of Electrical and Electronics Engineers (IEEE) originally for wireless transmission of high-definition video in the IEEE Std. 802.15.3c--2009 and then for wireless local area networks (WLANs) in the IEEE Std. 802.11ad--2012. Depending on the country-specific regulations, up to $6$ channels of $2.16$\,GHz each can be supported in IEEE 802.11ad terminals. The IEEE is currently working on further mmWave-specific features in IEEE P802.11ay, aiming to finalize it by the end of 2020.

The development of 5G-grade mobile mmWave communications is now almost complete, and vendors are testing their early implementations. The research focus is thus shifting towards determining the vision, requirements, and technology enablers for future wireless networks~\cite{beyond_5g_vision}. Due to the complexity in adopting mmWave in 5G, many stakeholders agree that the next-generation systems (often termed as 5G+) will primarily improve the efficiency of data exchange over existing $28$--$73$\,GHz bands. Simultaneously, harnessing the frequencies above $100$\,GHz is planned for 6G and beyond solutions as a key enabler for new applications demanding \emph{ultra-high user data rates exceeding tens of gigabits per second}~\cite{oulu_6g}.

The 6G vision is currently at the very early stages of its development. The support for futuristic scenarios is particularly considered, such as ubiquitous penetration of augmented and virtual reality (AR/VR) systems, holographic telepresence, and collective driving by autonomous robots~\cite{rangan_6G_use_cases_arxiv_commag}. From the technology perspective, the prospective 6G networks are also envisioned to introduce a number of novelties. These include full-duplex wireless communications, exploitation of machine learning for performance optimization, multi-radio multi-connectivity solutions, and many more~\cite{rap_6g}.

Concerning the extreme data rates to support, 6G systems are expected to complement microwave and mmWave connectivity options with wireless communications over the frequencies far above $100$\,GHz, particularly the Terahertz (THz) band communications ($300$\,GHz--$3$\,THz) and Visible Light Communications (VLC, $400$\,THz--$800$\,THz), both providing \emph{continuous bands of tens, hundreds, or even thousands of gigahertz}. The major players from academia, industry, standardization, and regulatory bodies have been actively exploring the feasibility of such communications over the last decades. For VLC, this effort has recently resulted in the creation of the IEEE Std. 802.15.7--2018 standard summarizing the physical (PHY) layer and the medium access control (MAC) sublayer for prospective optical wireless communications.

Standardization of THz wireless communications by the IEEE started in early 2008 when a ``Terahertz Interest Group'' (IGthz) has been established under the 802.15 umbrella. The key design choices and preliminary performance predictions made by the group members laid down the foundation of the IEEE Task Group on ``100G Wireless'' (TG100G, IEEE 802.15.3d) in 2014. The first IEEE standard for wireless communications over the low THz band (or \emph{sub-THz band}) by TG100G, IEEE Std. 802.15.3d--2017, has been officially approved in Fall 2017. While the key contributors are currently working on their prototype solutions by directly following the specifications~\cite{ieee_802_15_3d}, the successful adoption of the standard requires a better understanding of its main features and the existing limitations by a broader research community. The latter brings us to the main goal of this article.

In this paper, we present an overview of the IEEE Std. 802.15.3d--2017 -- the first IEEE-family standard for wireless communications over up to 69\,GHz-wide channels in the sub-THz band (specifically, $253$\,GHz--$322$\,GHz). We first review the target applications and usage scenarios in Section~\ref{sec:applications}. We then present the major details of the MAC and PHY layers in Section~\ref{sec:mac} and Section~\ref{sec:phy}, respectively. The initial performance assessment is reported in Section~\ref{sec:performance}. We finally outline the future work for the development and standardization of THz band communications in Section~\ref{sec:future}.

\section{IEEE~802.15.3d Motivation and Usage Scenarios}
\label{sec:applications}

The introduced IEEE Std. 802.15.3--2017 is an amendment to the IEEE Std. 802.15.3--2016. It defines a wireless PHY layer operating at $100$\,Gbit/s with a possible fallback to lower rates. The main objective is to illustrate the feasibility of sub-THz communications, as well as to provide a standard for low-complexity and high-rate connectivity. The rates are high enough for the emerging applications (see Fig.~\ref{fig:applications}), outlined further in this section in decreasing communication range.

The IEEE Std. 802.15.3d--2017, as one of the first standardization efforts for wireless communications at $300$\,GHz, primarily targets the applications that are feasible with the current level of electronics. Notably, the considered deployments are limited to point-to-point links between the devices that remain static or quasi-static during the data exchange~\cite{thz_ard}. The latter gives the key difference with the IEEE 802.11-family standards (particularly, IEEE Std. 802.11ad--2014 and IEEE P802.11ay) that primarily target mobile point-to-multipoint communications at $60$\,GHz with lower data rates.

The stationary and inherently point-to-point nature of the links allows reducing the complexity of the IEEE Std. 802.15.3d--2017 MAC since the requirements for multiple access and interference mitigation are relaxed. In addition, simpler procedures for initial access and device discovery can be used. Apart from the lower complexity in the implementation, the signaling overhead is reduced, allowing higher net data rates. Support of nodes mobility, as well as multiple channel access, is investigated for the future work~\cite{petrov_wcm_radar}.

\subsection{Wireless Fronthaul and Backhaul Links}
Future wireless networks are envisioned to utilize ultra-dense deployments of mmWave small cells to satisfy the projected capacity demands in crowded areas. One of the challenges here is to provide a reliable and high-rate \emph{backhaul} connection between the mmWave cell and the core network. In addition, when a small cell is equipped with several remote radio heads (RRHs), the stable and fast \emph{fronthaul} connection is needed between the RRHs and the base band unit (BBU).

Here, the massive utilization of optical fiber connections increases the ``transportation costs'' and complicates the deployment of mmWave small cells. Therefore, a wireless alternative has been proposed to complement the optical links in challenging setups. The approach is referred to as \emph{Integrated Access and Backhaul (IAB)} and suggests reusing the mmWave access technology for the backhaul~\cite{iab}. This solution is flexible but suffers from the limited data rate. The introduced IEEE Std. 802.15.3d--2017 combines the strengths of both options. On one side, the standard supports the data rates of $\approx100$\,Gbit/s -- that is an order of magnitude higher than the state-of-the-art mmWave technologies. On the other side, the use of IEEE Std. 802.15.3d--2017 notably simplifies the cabling.

\begin{figure}[!t]
 \centering
 \includegraphics[width=0.9\columnwidth]{pic/Fig1.pdf}
 \caption{Target applications for the IEEE 802.15.3d.}
 \label{fig:applications}
\end{figure}

\subsection{Additional Wireless Links in Data Centers}

Directional links via the IEEE Std. 802.15.3d--2017 can also be applied for wireless connections between the server racks in future data centers. One of the bottlenecks here is enabling efficient data exchange among the computing and storage nodes. Current solutions primarily utilize high-speed Ethernet via a fiber optic. However, the need for many cable channels challenges the reconfigurability and complicates cooling since some air streams become blocked.

Today, an alternative solution is actively studied, suggesting to utilize high-rate point-to-point wireless links when connecting the racks with each other. This approach enables a more flexible design of the data centers as well as reduces the amount of cabling. The IEEE Std. 802.15.3d--2017 aims to provide this ``inter-rack'' connectivity, thus complementing the fiber optics in certain deployments~\cite{thz_ard}.

\subsection{Kiosk Downloads}
Moving forward from business-oriented to consumer usage scenarios, the IEEE Std. 802.15.3d--2017 aims to offer high-rate data exchange for close-proximity communications. The two subcategories here are: (i) device-to-device (D2D) communications and (ii) kiosk download (see Fig.~\ref{fig:applications}).

The standardized sub-THz link offers an order of magnitude greater capacity in comparison with existing and emerging close-proximity wireless solutions in microwave or even mmWave bands. Particularly, the download time of a 2-hours-long movie ($900$\,MB) gets reduced from around $10$\,s with the IEEE 802.11ac to $0.1$\,s with the IEEE Std. 802.15.3d--2017~\cite{thz_ard}.

\subsection{Intra-Device Communications}
Further narrowing the range, not only consumer devices but also the electronic components inside a single device may be connected with the IEEE Std. 802.15.3d--2017. Today, modern computers are already equipped with several high-rate wired links interconnecting, e.g., the central processing unit (CPU) to the random-access memory (RAM), chipset to the network interface, and more. In addition, extremely-complex solutions are utilized in emerging massive multi-core CPUs to connect the computing cores and the shared cache memory.

The rapid growth in the data exchange inside the computers challenges their design. For instance, modern motherboards already have up to $12$ layers, while the emerging ``networks-on-chip'' for the future CPUs occupy over $30$\% of the processor space and power consumption. Here, the IEEE Std. 802.15.3d--2017 provides a possible alternative, enabling the high-rate data links between the critical components, while simultaneously simplifying the layout design~\cite{thz_ard,sergi_wnoc}.

\section{IEEE~802.15.3d MAC Layer}
\label{sec:mac}

In contrast to IEEE 802.11 and other standards for WLANs, IEEE Std. 802.15.3d--2017 support only point-to-point communications. The standard follows the MAC layer defined in the IEEE Std. 802.15.3e--2017~\cite{ieee_802_15_3e} and thus utilizes the concept of a \emph{pairnet} connecting no more than two devices. The point-to-point nature of the communications limits the range of possible use cases, but, at the same time, reduces the problem of interference and ``fighting for access''. The high-level description of the major signaling is presented in Fig.~\ref{fig:mac_signaling}.

\begin{figure}[!t]
  \centering
  \includegraphics[width=0.95\columnwidth]{pic/Fig2.pdf}
  \caption{Pairnet MAC layer signaling used in the IEEE 802.15.3d.}
  \label{fig:mac_signaling}
\end{figure}

The communications are divided into two major periods:
\begin{enumerate}
\item Pairnet setup period (PSP),
\item Pairnet associated period (PAP).
\end{enumerate}

During the PSP, one of the devices termed as the \emph{pairnet coordinator (PRC)} creates the pairnet and starts periodically sending \emph{Beacons} with the network details. The beacon frame includes the information on the number and duration of the access slots that can be used by the device willing to join the pairnet. Then, another device ready to joint the pairnet (termed as the \emph{pairnet device (PRDEV)}) responds to the beacon by sending the \emph{Association Request} at the beginning of one of the defined access slots. After a successful reception of the Association Request, the PRC stops transmitting beacons and sends an \emph{Association Response} to the PRDEV, ending the PSP. In addition to the pairnet setup, higher layer protocol setups, e.g., Internet Protocol (IP) or object exchange file transfer (OBEX), may be performed during the PSP. These additions are supported as both Association Request and Association Response frames have the fields to include higher-layer protocol information.

The second period -- PAP -- starts with the successful reception of the Association Response and is dedicated to the actual data exchange. During this period, both the PRC and the PRDEV transmit data frames and optional acknowledgments, starting with the PRDEV. The frames transmitted during PAP are separated by the defined \emph{short interframe space (SIFS)}. When one of the nodes wants to terminate the communications, it sends the \emph{Disassociation Request}. The PRC may also end PAP if there was no message received from the PRDEV during the defined timeout. In case the PRDEV does not want to terminate the communications but has no actual data, it may transmit a \emph{Probe Request} to restart the PRC timeout timer. Whenever the PAP is over (either via a Disassociation Request or by a timeout), the PRC switches back to the PSP: continues transmitting beacons and waiting for the new connections.

\section{IEEE~802.15.3d PHY Layer}
\label{sec:phy}

\subsection{General Overview}
\label{sec:phy_general}

The PHY layer in the IEEE Std. 802.15.3d--2017 has two modes:
\begin{itemize}
\item THz single carrier mode (\emph{THz-SC PHY}), as described in Subsection~\ref{sec:phy_sc},
\item THz on-off keying mode (\emph{THz-OOK PHY}), as described in Subsection~\ref{sec:phy_ook}.
\end{itemize}

The first of these PHY modes, THz-SC PHY, is designed for the high data rates. This mode primarily targets the bandwidth-oriented use cases from Section~\ref{sec:applications}, such as wireless fronthaul/backhaul and additional links in the data center. On its turn, the second mode, THz-OOK PHY, targets the lower-cost sub-THz devices that are not able to utilize complex signals and thus have to rely on the amplitude information. This is the case, for example, with solutions employing resonant tunneling diodes (RTD). Still, the rates of tens of Gbit/s can be achieved with the THz-OOK PHY when the widest channels are used.

\subsection{Supported Channels}
\label{sec:channels}
The standard may operate over the sub-THz frequencies between $252.72$\,GHz and $321.84$\,GHz. In total, there are $69$ overlapping channels (see Fig.~\ref{fig:channels}). There are $8$ supported channel bandwidths ranging from $2.16$\,GHz, as in mmWave IEEE Std. 802.11ad--2012, and up to massive $69$\,GHz. Depending on the use case and the hardware capabilities, either the entire frequency range can be allocated for a single $69.12$\,GHz-wide channel (CHNL\_ID=69) or shared between several smaller channels. The bandwidths of the smaller channels are integer multiples of $2.16$ GHz. The channel number $41$ with the bandwidth of $4.32$\,GHz ($2\times2.16$\,GHz) is the default channel.

Following the decision at World Radio Conference 2019 (WRC-2019), all the bands depicted in Fig.~\ref{fig:channels} are available for THz communications globally if specific conditions to protect radio astronomy and earth exploration satellite service are met, as further discussed in Section~\ref{subsec:regulatory}. These conditions do not explicitly specify any transmit power limits and are applicable, in practice, primarily to the narrow areas surrounding the ground radio astronomy stations.

\begin{figure}[!h]
 \centering
 \includegraphics[width=1.0\columnwidth]{pic/Fig3.pdf}
 \caption{IEEE 802.15.3d channel plan.}
 \label{fig:channels}
\end{figure}

\begin{figure}[!t]
 \centering
 \includegraphics[width=1.0\columnwidth]{pic/Fig4.pdf}
 \caption{IEEE 802.15.3d frame format.}
 \label{fig:frame_header}
\end{figure}

\subsection{THz-SC PHY Layer Mode}
\label{sec:phy_sc}
The THz-SC mode is primarily intended for high-rate communications. According to~\cite{ieee_802_15_3d}, the frame length varies from $2048$ bytes to $2$\,$099$\,$200$ bytes, not including the PHY preamble or the base header. THz-SC PHY supports $6$ different modulations. The first four are phase shift keying modulations: binary (BPSK), quadrature (QPSK), 8-Phase (8-PSK), and 8-Phase Asymmetric (8-APSK). The mode also supports quadrature amplitude modulation with $16$ and $64$ constellation points: 16-QAM and 64-QAM. Forward error correction (FEC) is performed via one of the two low-density parity-check (LDPC) codes: either a high-rate 14/15 LDPC (1440,1344) or a low-rate 11/15 LDPC (1440,1056). BPSK and QPSK modulations are mandatory for the THz-SC mode, the others are optional.

The frame format used in the THz-SC mode is illustrated in Fig.~\ref{fig:frame_header}. The frame starts with a PHY preamble to aid frame detection, timing acquisition, and channel estimation at the receiver. There are two options for the PHY preamble: a long version is transmitted during the PSP, while the short one is transmitted during the PAP. The preamble is followed by a PHY header, detailed in Fig.~\ref{fig:frame_header}. Both the MAC and the PHY headers are protected with the header check sequence (HCS), implemented using the cyclic redundancy check (CRC), particularly, CRC-16. To increase the robustness, the combination of the PHY header, MAC header, and HCS is encoded to concatenated codewords of an extended Hamming (EH) code. Additional header stuff bits may be added to get an integer length of the data portion in the block. Finally, the optional header stuff bits are followed by the frame payload.

\subsection{THz-OOK PHY Layer Mode}
\label{sec:phy_ook}
In contrast to the THz-SC mode aiming for the highest possible data rates, the THz-OOK mode of the IEEE~802.15.3d PHY layer is intended for simpler devices. The primary use cases for this mode are kiosk downloading and intra-device communications, as described in Section~\ref{sec:applications}. The target transmission range for this mode thus does not exceed a few meters. The THz-OOK mode supports a single low-complexity modulation scheme -- on-off keying (OOK) -- where the data are represented either as the presence (logical ``1'') or absence (logical ``0'') of a carrier wave. On top of the OOK modulation, three FEC schemes are supported. The (240,224)-Reed Solomon (RS) code is mandatory for this mode, as it allows simple decoding without soft decision information. The other two FEC options are the LDPC-based schemes outlined in the previous subsection. They are \emph{optional} for the THz-OOK mode and enable the use of soft-decision information.

The frame structure used in THz-OOK mode is similar to the one used in THz-SC mode (see Fig.~\ref{fig:frame_header}) with just the minor differences. Particularly, as there is just a single modulation (OOK) plus three FEC options, only two bits are used in the PHY header to define the current modulation and coding scheme (MCS), instead of four bits in THz-SC. Other fields are identical, thus making the THz-OOK PHY header two bits shorter than the one for the THz-SC mode. The target minimal receiver sensitivity levels for these two modes are also the same and equal to $-67$\,dBm (with 11/15 LDPC code and the smallest bandwidth of $2.16$\,GHz).

\section{Initial Performance Evaluation}
\label{sec:performance}

\begin{table}[!b]
\caption{Key numerical study parameters}
\label{tab:params}
\begin{center}
\begin{tabular}{|p{0.24\columnwidth}|p{0.135\columnwidth}|p{0.125\columnwidth}|p{0.135\columnwidth}|p{0.125\columnwidth}|}
\hline
 & \textbf{Fronthaul /backhaul} & \textbf{Data \newline center} & \textbf{Kiosk download} & \textbf{Intra-device}\\
\hline
Tx power & 25 dBm & 10 dBm & 0 dBm & 0 dBm\\
\hline
Tx gain & 30 dB & 30 dB & 24 dB & 6 dB\\
\hline
Rx gain & 30 dB & 30 dB & 12 dB & 6 dB\\
\hline
\hline
Carrier frequency & \multicolumn{4}{c|}{300 GHz}\\
\hline
Rx noise figure & \multicolumn{4}{c|}{8 dB}\\
\hline
Error Vector & \multicolumn{4}{l|}{$[-22 ... -3]$ dB, depending on the selected MCS}\\
Magnitude (EVM) & \multicolumn{4}{l|}{and channel bandwidth as detailed in~\cite{thz_simulator}}\\
\hline
\end{tabular}
\end{center}
\end{table}

In this section, we report the results of the initial performance evaluation of the IEEE Std. 802.15.3d--2017.

\textcolor{black}{We first recall the maximum PHY-layer data rates reported for the IEEE Std. 802.15.3d--2017 (particularly, Table 13-4 for the THz-SC mode and Table~13-12 for the THz-OOK mode~\cite{ieee_802_15_3d}).} Point-to-point connectivity is assumed with the SNR sufficiently high to guarantee a negligible bit error rate (BER) with the selected MCS. Following~\cite{ieee_802_15_3d}, the data rates with the THz-OOK mode vary from $1.64$\,Gbit/s for the smallest band of $2.16$\,GHz up to $52.56$\,Gbit/s for the largest $69.12$\,GHz-wide channel. The more advanced THz-SC can theoretically support higher rates up to $315.39$\,Gbit/s with 64-QAM modulation and\,14/15\,LDPC\,code.

We then proceed with a study on the communication range, as reported in Fig.~\ref{fig:comm_range}. This figure illustrates the maximum distance between the PRC and PRDEV that still allows for reliable data exchange as a function of the channel bandwidth. For this study, a link-level simulator was developed, as detailed in~\cite{thz_simulator}. For each of the use cases, the simulator utilizes the corresponding application-specific parameters and channel models for sub-THz communication introduced in~\cite{thz_cmd}. The efficient 64-QAM modulation is used for the long-range setups (fronthaul/backhaul and data center). A more robust 8-APSK modulation is utilized for kiosk downloading and intra-device. For all the use cases, the high-rate 14/15 LDPC code is applied. The key modeling parameters are given in Table~\ref{tab:params}.

Following Fig.~\ref{fig:comm_range}, the range vary from centimeters for the intra-device and up to several hundred meters for the fronthaul/backhaul. At the same time, as detailed in~\cite{thz_simulator}, the use of highly-directional antennas is mandatory on both PRC and PRDEV sides, especially for the long-range fronthauling and backhauling. The antenna gains of $\geq30$\,dB are required for this use case to mitigate the spatial loss and the absorption. Fig.~\ref{fig:comm_range} also shows the maximum range, while still maintaining a $100+$\,Gbit/s data rate. These values reach $\approx100$\,m for the fronthaul/backhaul, $17$\,m for the data center, $0.61$\,m for the kiosk downloading, and $0.03$\,m ($3$\,cm) for the intra-device connectivity. The results are several times lower than the values achievable with a $2.16$\,GHz channel, highlighting a trade-off between the range and the link-level performance. In practical implementations, operators may adjust this range--rate balance depending on the particular application and deployment.

Summarizing, the preliminary performance assessment of the IEEE Std. 802.15.3d--2017 illustrates that the essential rate- and range-related requirements are satisfied in all the described use cases. Nevertheless, more detailed evaluations are still in progress, aiming to refine the results in specific practical setups, such as wireless backhauling in realistic environments and the data centers with complex topology~\cite{bh_detailed}.

\begin{figure}[!t]
 \centering
 \includegraphics[width=1.0\columnwidth]{pic/Fig5.pdf}
 \caption{Range and rate values achievable in target IEEE 802.15.3d use cases.}
 \label{fig:comm_range}
\end{figure}

\section{The Road Ahead in Standardizing\\THz Band Communications}
\label{sec:future}

\subsection{Regulatory Aspects of THz Band Communications in 6G}
\label{subsec:regulatory}
THz communications is one of the candidates to enable ultra-high data rates in 6G networks~\cite{rap_6g}. THz systems cannot maintain wide coverage, and therefore they will not suit as the general-purpose 6G air interface. Still, THz connections will be an important component providing extremely-high-rate data pipes. To facilitate this, an unprecedented amount of $160$\, GHz of the spectrum was recently opened for THz communications after the WRC--2019, which has included a new note RR No. FN5.564A describing the conditions to use the $275$\,GHz--$450$\,GHz band by land mobile and fixed service~\cite{wrc_final}.

One of the major issues analyzed in the regulatory process is the potential interference from the THz communications to existing ground-to-orbit links. The effect has been extensively studied over recent years. As a result, sharing with passive services is still necessary within certain parts of the specified band. In contrast, within $252$--$296$\,GHz, $306$--$313$\,GHz, $318$--$333$\,GHz, and $356$--$450$\,GHz no conditions are necessary to protect earth exploration-satellite service (EESS), and well manageable conditions exist to protect the radio astronomy service (RAS)~\cite{wrc_final}. Hence, massive bands are already available for future THz communication systems.

\textcolor{black}{The presented WRC allocations are global, they hold for the entire world and provide the baseline for further regulatory activities. No specific limits for the transmit power levels are specified to date. Additional region-specific and national regulations are expected to appear in the coming years. These regulations will follow the general agreements summarized in~\cite{wrc_final} but might specify more details, e.g., determine the allowed power levels, the allocations of the sub-bands, and the necessity and nature of the required coexistence mechanisms.
}

\subsection{Prototyping Activities for IEEE Std. 802.15.3d}
The implementation of the end-to-end communication system operating in the IEEE Std. 802.15.3d--2017 bands is currently an active work in progress. These demonstrators are primarily built within large-scale research and infrastructure projects. Particularly, the demonstrator developed within the EU-Japan project THoR applies the chipsets following the implementation in~\cite{ieee_802_15_3e}, which uses the same MCSs as defined in the THz-SC PHY mode. Another solution based on RTDs has been presented in the EU H2020 project iBROW, aiming to provide low-complexity and cost-efficient implementation that is more suitable for the THz-OOK mode.

Started several years ago, these prototyping activities have been pushed further by the final approval of the IEEE Std. 802.15.3d--2017 specification, as well as the recent frequency allocations made at the WRC-19~\cite{wrc_final}. A recent example here is a transceiver chip capable to utilize the channels 49--51 and 66 from Fig.~\ref{fig:channels}~\cite{extra_cmos_thz_tranceiver}. Further prototyping activities are upcoming, following the announced decisions and the availability of the hardware needed for sub-THz radio links.

\subsection{Towards Prospective THz Band WLANs}
The IEEE Std. 802.15.3d--2017 specifies wireless THz connections for fixed point-to-point links only, where the (i) directions of the antennas at both ends of the links are known; (ii) interference can be mitigated by proper links planning~\cite{bh_detailed}; and (iii) no fight for access exists. In prospective THz WLANs, these conditions do not hold. Therefore, the proper solutions have to be developed first to address the listed aspects.

Once scientific breakthroughs in these areas are made, the corresponding algorithms have to be included in standards. Several possibilities are currently under consideration in the IEEE 802.15 TAG THz and other teams, such as developing amendments to the existing IEEE Std. 802.15.3d--2017, incorporating selected 802.15.3d features into future releases of IEEE 802.11, or creating completely new standards. Today, the IEEE 802.15 TAG THz is also exploring novel solutions for physical and link layers to operate over $10$--$100$\,GHz-wide channels in sub-THz and THz bands. On this way forward, the ideas and solutions delivered in the IEEE Std. 802.15.3d--2017 will contribute to the appearance of THz communications, as well as form the basis for prospective 6G wireless systems.

\balance


\section{Introduction}

High transition temperature ($T_c$) cuprate superconductors remain a central topic of condensed matter research since their discovery in the 1980's.
Among their many fascinating and mysterious properties\cite{Proust}, the existence of pseudogaps and Fermi arcs in the underdoped regions of their phase diagrams is particularly intriguing\cite{Yoshida}. Both of them indicate the breakdown of the standard Fermi liquid\cite{Book} description of the non-superconducting state above $T_c$, on which the Bardeen-Cooper-Schrieffer (BCS) theory\cite{Book} of the superconducting state is founded. Such non-Fermi liquid behavior is also seen in many other properties including lack of coherence in quasiparticle excitations and unusual transport properties in the normal state\cite{Proust,Yoshida}, which are likely related to the physics of neighboring Mott insulator phase. Undertanding such non-Fermi liquid physics, including (and in particular) the origin of Fermi arcs and pseudogap, is an exciting challenge we face.

In this paper we present an extremely simple and exactly solvable model, and show that Fermi arcs and pseudogap appear very naturally (and hand-in-hand) in its ground state and excitation spectrum. Our model is a variant of a model introduced by Hatsugai and Kohmoto (HK)\cite{hk} (a model similar to that of HK was considered earlier by Baskaran\cite{Baskaran}).
We note the HK model is also the basis of a recent theory of superconductivity in doped Mott insulators\cite{Phillips}, building on earlier phenomenological work\cite{Setty}. We would like to emphasize upfront that it is {\em not} our purpose to present a comprehensive theory for the cuprate phenomenology based on such a simple model, which, as we will explain in more detail below, misses certain important aspects of the cuprate physics. Its value is providing a proof of principle, namely Fermi arcs and pseudogap can exist. Given their importance, we believe a simple model that provides an existence proof is highly desirable. It may also provide a starting point for building more realistic models that can eventually lead to a comprehensive understanding of the extremely rich phenomenology of the cuprates.

\section{Model and its solution}

We consider the following model:
\begin{equation}
H=\sum_{{\bf k}}[\epsilon_{{\bf k}}(n_{{\bf k}\uparrow} + n_{{\bf k}\downarrow})+u_{{\bf k}} n_{{\bf k}\uparrow}n_{{\bf k}\downarrow}],
\label{eq:model}
\end{equation}
where the single particle energy $\epsilon_{{\bf k}}$ is measured from the chemical potential, and $n_{{\bf k}\sigma}$ is the fermion occupation for momentum ${\bf k}$ and spin $\sigma = \uparrow,\downarrow$. For simplicity and concreteness, we consider, {\em e.g.}, a square lattice with (anisotropic) nearest neighbor hopping, resulting in
\begin{equation}
\epsilon_{{\bf k}} =-2t_x\cos(k_x) - 2t_y\cos(k_y) -\mu,
\label{eq:dispersion}
\end{equation}
where $t_x$ and $t_y$ are hopping along x- and y-directions, and $\mu$ is chemical potential.
If the interaction energy $u_{{\bf k}} = U$ is a constant, (\ref{eq:model}) reduces to the HK model\cite{hk}. The momentum dependence of $u_{{\bf k}}$ breaks the very special center of mass conservation of the HK interaction, and can be understood as renormalization of dispersion of the interacting pair. Since such interaction is expected to emerge at low energy from the interplay of single particle dispersion and more realistic two-body interaction, existence of such dispersion renormalization is physically very natural. Again for simplicity and concreteness, we consider
\begin{equation}
u_{{\bf k}} =-2T_x\cos(k_x) - 2T_y\cos(k_y) +U,
\label{eq:pairdispersion}
\end{equation}
where $T_x$ and $T_y$ may be understood as the renormalization of hopping of center of mass of the interacting pair.
The solvability of the model lies in the fact that $n_{{\bf k}\uparrow}$ and $n_{{\bf k}\downarrow}$ are conserved quantities. We emphasize the details of the momentum dependence of Eqs. (\ref{eq:dispersion}) and (\ref{eq:pairdispersion}) are {\em not} important for our discussions below.

If $u_{{\bf k}}$ is either positive or negative definite, the solution is qualitatively very similar to that of the HK model. In the following we show the situation is very different when $u_{{\bf k}}$ changes sign in momentum space. In particular when the $u_{{\bf k}}=0$ surface intersects with the non-interacting Fermi surface $\epsilon_{{\bf k}} = 0$, the ground state supports Fermi arcs.

Defining $n_{{\bf k}} = n_{{\bf k}\uparrow} + n_{{\bf k}\downarrow}$, we find in the ground state
\begin{eqnarray}
n_{{\bf k}}= \left\lbrace\begin{array}{cc}
0, &\epsilon_{{\bf k}} > 0 \hskip 0.2 cm {\rm and} \hskip 0.2 cm 2\epsilon_{{\bf k}} + u_{{\bf k}} > 0;\\
1, &\epsilon_{{\bf k}} < 0 \hskip 0.2 cm {\rm and} \hskip 0.2 cm \epsilon_{{\bf k}} + u_{{\bf k}} > 0;\\
2, &\epsilon_{{\bf k}} + u_{{\bf k}} < 0 \hskip 0.2 cm {\rm and} \hskip 0.2 cm 2\epsilon_{{\bf k}} + u_{{\bf k}} < 0. \end{array}\right.
\label{eq:occupation}
\end{eqnarray}
The surfaces (or lines) in momentum space across which $n$ jumps are given by
\begin{eqnarray}
\label{eq:boundary 1}
&&\epsilon_{{\bf k}} = 0,\\
\label{eq:boundary 2}
&&\epsilon_{{\bf k}} + u_{{\bf k}} = 0,\\
\label{eq:boundary 3}
&& 2\epsilon_{{\bf k}} + u_{{\bf k}} = 0.
\label{eq:boundaries}
\end{eqnarray}
Obviously all three lines intersect with the loop
\begin{eqnarray}
u_{{\bf k}}=0
\label{eq:boundary 4}
\end{eqnarray}
at the same points. In Fig. 1(a) all four lines are plotted in a simple representative case, resulting in the occupation pattern Fig. 1(b).

\begin{figure
\center
\includegraphics[width=0.45\textwidth]{fig1.pdf}
\caption{(a) Solutions of Eqs. (\ref{eq:boundary 1}) (black line), (\ref{eq:boundary 2}) (red line), (\ref{eq:boundary 3}) (blue line) and (\ref{eq:boundary 4}) (green line). We assume the quantities that vanish along these loops are negative inside and positive outside. (b) Occupation pattern corresponding to panel (a). The black and red lines are Fermi arcs while the blue lines are pseudo Fermi surfaces along which there is a pseudogap.}
\label{fig:fig1}
\end{figure}

We now distinguish two different types of boundaries in Fig. 1(b).

$\bullet$ $\Delta n_{{\bf k}} =2$. This case looks (superficially) like the Fermi surface of the non-interacting Fermi gas, but is {\em not} a Fermi surface, because single fermion excitations are actually {\em gapped} there with gap $|u_{{\bf k}}|/2$.
As a result we refer to such boundaries pseudo Fermi surfaces, and identify the fermion gap there the pseudogap (to be distinguished from superconducting gap; more on this point later).

$\bullet$ $\Delta n_{{\bf k}} =1$. This case corresponds to our Fermi arcs, and are Fermi surfaces in the sense single fermion excitations are gapless there.

Having identified the Fermi arcs and pseudogap regions in our model, a few comments are in order.

$\bullet$
As pointed out by HK\cite{hk} the region with $n_{{\bf k}}=1$ (surround by the Fermi arcs in our case) has a massive degeneracy associated with a two-fold degeneracy for each ${\bf k}$. This is an artifact of the model. In fact this degeneracy can be easily removed by introducing a Zeeman splitting $\Delta_Z$ between up- and down-spin fermions:
\begin{eqnarray}
H_Z&=&\Delta_Z\sum_{{\bf k}}n_{{\bf k}\uparrow}\\
&=&\frac{\Delta_Z}{2}\sum_{{\bf k}}(n_{{\bf k}\uparrow} - n_{{\bf k}\downarrow})+\frac{\Delta_Z}{2}\sum_{{\bf k}}(n_{{\bf k}\uparrow} + n_{{\bf k}\downarrow}),
\label{eq:Zeeman}
\end{eqnarray}
and our previous analysis goes through exactly the same way, as long as one replaces $u_{{\bf k}}$ by $u_{{\bf k}} + \Delta_Z$ (note the last term above results in a shift of chemical potential $\mu$). Then the two Fermi arcs surrounding the $n_{{\bf k}}=1$ region are Fermi surfaces for the up- and down-spin fermions respectively. Even if $\Delta_Z=0$, as we will demonstrated later, ferromagnetic (Stoner) instability will remove this degeneracy when generic weak repulsive interaction is present. For the ease of discussion we will assume such a splitting is present (regardless of its origin), although most of our conclusions hold even when the massive degeneracy mentioned above is present.

$\bullet$ The pseudogap is clearly a singlet pairing gap. However it is {\em not} due to a fluctuating BCS pairing order parameter. In fact to (literally) the opposite, the pairing here is between up- and down-spin fermions at the {\em same} momentum, instead of opposite momenta as in BCS. The pairing here could perhaps be viewed as an extreme version of fluctuating pair density wave (PDW)\cite{Agterberg}, in which {\em every} pair has a different (and very large) momentum. We note a fluctuating PDW has been proposed by Lee and co-workers\cite{Lee} to be the mother of many competing states in underdoped cuprates, and hold the key to Fermi arcs and pseudogap.

$\bullet$ As discussed earlier, Zeeman splitting results in an effective increase of $u_{{\bf k}}$. This widens the Fermi arcs at the expense of the pseudo Fermi surface (or pseudogap region). More quantitatively, the pseudogap as well as pseudo Fermi surface are eliminated when
\begin{eqnarray}
\Delta_Z = |u_{\rm min}|=\Delta_{\rm pseudo},
\label{eq:Zeeman}
\end{eqnarray}
where $u_{\rm min}$ is the (negative) minimum value of $u_{{\bf k}}$ along the pseudo Fermi surface, whose magnitude is the maximum value of pseudogap $\Delta_{\rm pseudo}$.
Such suppression of pseudogap by Zeeman splitting, and in particular the simple linear scaling with $\Delta_Z$, have indeed been seen in the cuprates\cite{Shibauchi}.

\section{Angle-resolved photoemission spectroscopy and Luttinger sum rule}

An angle-resolved photoemission spectroscopy (ARPES) experiment\cite{Yoshida} measures the negative energy portion ($\omega < 0$, representing states occupied by electrons) of the electron spectral function, which is the imaginary part of the retarded electron Green's function:
\begin{eqnarray}
G(\omega, {\bf k}) = \frac{1}{\omega - \epsilon_{\bf k} + i\delta}
\end{eqnarray}
if $n_{\bf k} = 0$ in the ground state,
\begin{eqnarray}
G(\omega, {\bf k}) = \frac{1}{\omega - \epsilon_{\bf k} - u_{\bf k} + i\delta}
\end{eqnarray}
if $n_{\bf k} = 2$ in the ground state, and
\begin{eqnarray}
G_\uparrow(\omega, {\bf k}) &=& \frac{1}{\omega - \epsilon_{\bf k} - u_{\bf k} + i\delta},\\
G_\downarrow(\omega, {\bf k}) &=& \frac{1}{\omega - \epsilon_{\bf k} + i\delta}
\end{eqnarray}
if $n_{\bf k} = 1$ in the ground state. In the above $\delta$ is an infinitesimal positive, and we took the limit $\Delta_Z\rightarrow 0$ so that we do not need to distinguish between up- and down-spin Green's function in regions with $n_{\bf k} = 0$ and $n_{\bf k} = 2$ in the ground state. None of the conclusions below depends on this simplification. Due to the single pole structure, the spectral function contains a single delta-function peak in our model.

ARPES will find the Fermi arcs exactly the same way ordinary Fermi surfaces are identified. It will (in principle) see both branches (and therefore a closed surface) if it is not spin-resolved, although not necessarily resolving them if their separation is smaller than the momentum resolution of ARPES. On the other hand a spin-resolved ARPES will see one branch but not the other, resulting in open arcs.

On the other hand ARPES will find a gap along the pseudo Fermi surface, which, for $|u_{{\bf k}}|$ small compared to Fermi energy, will be very close to the Fermi surface of non-interacting fermions. ARPES cannot distinguish this gap from the superconducting gap; since there is no superconducting order in our ground state, it should be identified as a pseudogap. Indeed it suppresses the density of states and various thermodynamic response functions compared with non-interacting fermions or a Fermi liquid.

It should be emphasized that ARPES finds very broad electron spectral functions in the underdoped cuprates, with very small weight in a coherent quasiparticle peak (if present)\cite{Yoshida}. This means there is very strong scattering that render single particle excitations incoherent. This important piece of physics, which is related to the fact that underdoped cuprates are doped Mott insulators, is missing in our model.

Since our Fermi arcs are {\em not} closed, Luttinger's theorem (also known as Luttinger sum rule)\cite{AGD} is not satisfied the way in a Fermi liquid, where the Fermi surfaces are closed. This is, of course, just one of many indications that we have a non-Fermi liquid here. The way Luttinger sum rule gets satisfied here involves the combination of the Fermi arcs {\em and} the pseudo Fermi surfaces. This is because the fermion Green's function at zero energy, $G(\omega=0, {\bf k})$ changes sign across the pseudo Fermi surfaces (just like real Fermi surfaces), even though the fermion spectral weight there is {\em zero} (instead of a delta-function peak as for real Fermi surfaces of a Fermi liquid, and reflecting the pseudogap). Reversing the logic, Luttinger sum rule {\em requires} the co-existence of open Fermi arcs and pseudo Fermi surface(s) of the type discussed here.

\section{Quantum Oscillations}

The conclusions above follow straightforwardly from the exact solution of (\ref{eq:model}). In this section we discuss a more subtle question, namely how the system responds to a perpendicular magnetic field, which gives rise to quantum oscillations and represents another important probe of the pseudogap and Fermi arcs\cite{Proust,Sebastian}. Here we focus on the orbital coupling of the magnetic field, and ignore its (experimentally much weaker) Zeeman effect which, as discussed earlier, is straightforward to include.

We have three types of low-energy excitations in our model: (i) Spin-up fermions, whose dispersion is
\begin{eqnarray}
E^\uparrow_{{\bf k}}=\epsilon_{{\bf k}},
\end{eqnarray}
living near the black Fermi arcs of Fig. 1(b); (ii) Spin-down fermions, whose dispersion is
\begin{eqnarray}
E^\downarrow_{{\bf k}}=\epsilon_{{\bf k}} + u_{{\bf k}},
\end{eqnarray}
living near the red Fermi arcs of Fig. 1(b); and (iii) Singlet fermion pairs, whose dispersion is
\begin{eqnarray}
E^{\rm pair}_{2{\bf k}}=2\epsilon_{{\bf k}} + u_{{\bf k}},
\end{eqnarray}
living near the pseudo Fermi surfaces (blue) of Fig. 1(b). In the presence of a magnetic field ${\bf B}$, each of them propagate around their respective iso-energy contour following the semi-classical equations of motion\cite{Book}:
\begin{eqnarray}
\frac{d{\bf k}}{dt} &=& -\frac{q}{\hbar c}{\bf v}_{{\bf k}}\times {\bf B},\\
\hbar {\bf v}_{{\bf k}} &=& \nabla_{{\bf k}}(E_{{\bf k}}),
\label{eq:oscillation}
\end{eqnarray}
where $q=e$ for single fermions and $q=2e$ for the singlet pair.
If there were no interaction/restriction among them, these would be the contours illustrated in Fig. 1(a). Quantizing such periodic motion gives rise to Landau levels and density of states oscillating periodically with $1/B$ with a periodicity determined by the area of these contours\cite{Book}. The situation here is quite different: As each of them arrives at the point where $u_{{\bf k}}=0$, continuing with the trajectory given by their dispersion would take them into a region forbidden either by energy conservation or Pauli blocking. This is quite similar to a quantum mechanical particle following a semi-classical trajectory hitting an energy barrier. The resolution is, of course, that the particle quickly tunnels through the classically forbidden region and lands at the closest classically allowed point (due to the chirality imposed by the magnetic field back-scattering is also forbidden), and then continues with the semi-classical motion. This results in the momentum space trajectories illustrated in Fig. 2. Since the tunneling process is expected to be very quick, the periods of such motions are much shorter because the particles traverse only a small part of their semi-classical trajectories. As a result the quantum oscillations resulting from such trajectories have much smaller effective ``areas". If interpreted in the same way as quantum oscillations in Fermi liquids, one would conclude there exist Fermi ``pockets" that are much smaller than the original Fermi surface satisfying the Luttinger sum rule. This is indeed what is seen in the pseudogap region of the cuprate phase diagram\cite{Sebastian}.

\begin{figure
\center
\includegraphics[width=0.4\textwidth]{fig2.pdf}
\caption{The momentum space trajectories of the three types of low-energy excitations of our model. Solid lines represent semi-classical trajectories while dotted lines represent tunneling through forbidden region.}
\label{fig:fig2}
\end{figure}

We note among the many competing theories of small Fermi pockets\cite{Proust}, ours belongs to the group in which there is {\em no} broken translation symmetry (and the corresponding Fermi surface folding). Within this group the theory of Ref. \cite{Pereg-Barnea} is perhaps the closest to ours in spirit. We would like to emphasize Ref. \cite{Pereg-Barnea} models the pseudogap by a superconducting gap (which is argued to lose phase coherence at long length scales so that the system is {\em not} a superconductor); in our model the pseudogap has nothing to do with superconducting gap, and (as noted earlier) can be viewed as due to pairing fluctuation of the PDW type instead.

\section{Robustness of results and possible superconducting instability}

In this section we demonstrate the robustness of the results that follow from the exact solution of (\ref{eq:model}) against more generic repulsive interactions. To this end we consider a perturbing Hubbard interaction:
\begin{eqnarray}
H' = J \sum_{i}n_{i\uparrow}n_{i\downarrow} = \frac{J}{N}\sum_{{\bf k}{\bf k}'{\bf q}}c^\dagger_{{\bf k}+{\bf q}\uparrow}c^\dagger_{{\bf k}'-{\bf q}\downarrow}c_{{\bf k}'\downarrow}c_{{\bf k}\uparrow},
\label{eq:Hubbard}
\end{eqnarray}
where $i$ is a site index, $J$ is Hubbard interaction strength (the unusual notation is going to be justified soon) and $N$ is system size. To linear order in $J$, we only need to evaluate $H'$ in the Fock subspace in which $H$ of (\ref{eq:model}) is degenerate:
\begin{eqnarray}
H' = (J/N)\sum_{({\bf k}{\bf k}')}[f(n_{\bf k}, n_{{\bf k}'}) - 2 {\bf S}_{\bf k}\cdot {\bf S}_{{\bf k}'}],
\label{eq:1st order Hubbard}
\end{eqnarray}
where ${\bf S}_{\bf k}$ is the spin operator associated with momentum ${\bf k}$, taking values $s=1/2$ for $n_{\bf k} =1$ and $s=0$ for $n_{\bf k} =0, 2$; $f(0, 0)=f(0, 1)=f(1,0)=0$, $f(1,1)=1/2$, $f(1,2)=f(2,1)=1$, and $f(2,2)=2$. The summation is over all pairs $({\bf k}{\bf k}')$. The ferromagnetic coupling in (\ref{eq:1st order Hubbard}) for $J > 0$ results from the {\em direct} exchange in momentum space; this is opposite to what happens in the usual Hubbard model at half-filling, where a super exchange in real space results in an antiferromagnetic coupling\cite{Book}. As anticipated earlier this ferromagnetic coupling fully polarizes the fermion spins in the $n_{\bf k} = 1$ region and removes the massive ground state degeneracy of (\ref{eq:model}). Due to the long-range nature of the spin-spin interaction (in momentum space) in (\ref{eq:1st order Hubbard}), there are {\em no} gapless spin-wave excitations here. The origin of such unusual behavior can be traced to the long-range nature of the $u_{\bf k}$ interaction (in {\em real} space) in (\ref{eq:model}). As a result this spontaneous magnetization has essentially the same effect as a Zeeman splitting. We thus expect our earlier results are robust for $J > 0$.

The situation is quite different for $J < 0$ (attractive Hubbard interaction). In this case the spin-spin interaction in (\ref{eq:1st order Hubbard}) is antiferromagnetic, and highly frustrated due to its long-range nature. Due to this frustration any state with ${\bf S}_{\rm tot} = \sum_{\bf k}{\bf S}_{\bf k} = 0$ is a ground state, and the massive ground state degeneracy is reduced but {\em not} lifted at linear order in $J$, and remains exponentially large in system size. The fate of the ground state is thus determined by higher-order processes in $J$. To this end we expect the BCS pairing interaction considered in Ref. \cite{Phillips} to dominate, resulting in a superconducting ground state, in whose spectrum the superconducting gap coexists with the pseudogap discussed earlier. We leave a detailed analysis of this case to future work.

\section{Summary and Discussions}

To summarize, we have presented a simple model that gives rise to Fermi arcs and pseudogap, in ways that are qualitatively consistent with phenomenology of underdoped cuprates, especially those of ARPES and quantum oscillation experiments. As emphasized earlier it is not our intention to present a comprehensive theory of cuprate phenomenology here. On the other hand we believe such a simple model is interesting in its own right, and may well contain the basic ingredients of the strong correlation physics responsible for the relevant phenomenology of cuprates and other systems.

We could attempt a more quantitative comparison with cuprate phenomenology, by choosing appropriate momentum, doping and perhaps even temperature dependence of parameters of our model. For example Fig. 1 corresponds to $\epsilon_{\bf k}$ and $u_{\bf k}$ with a $C_2$ symmetry, instead of $C_4$ symmetry that is appropriate to (some of the) cuprates. It is straightforward to enlarge their symmetry, which will result in four instead of two pairs of Fermi arcs. One can in principle force very detailed agreement, by fine-tuning $\epsilon_{\bf k}$ and $u_{\bf k}$ in ad hoc ways. On the other hand as stated before the main point of the present work is providing a proof of principle, and possible starting point for future work. Fine tuning would not add any additional insight, but might dilute our main message instead. On that note let us reiterate that there is also important cuprate physics missing in our model, especially those associated with the Mott insulator. In addition to those already mentioned we add that antiferromagnetism plays a very important role in cuprates, while in our model it appears to be suppressed; instead ferromagnetism is found with generic repulsive interactions. Understanding the possible interplay between Fermi arcs/pseudogap and Mott physics in the context of our model is a very interesting direction for future work.

\section*{Acknowledgments}
The author thanks Chandan Setty for correspondences that motivated the present work, and Joseph Ghobrial for graphics assistance.
This work was supported by the National Science Foundation Grant No. DMR-1932796, and performed at the National High Magnetic Field Laboratory, which is supported by National Science Foundation Cooperative Agreement No. DMR-1644779, and the State of Florida.


\section{Introduction}
Finding roots of univariate complex polynomials numerically is one of the fundamental problems for numerical
algebraic geometry. %
Many numerical algorithms for finding roots can be interpreted as complex dynamical systems.
For example, the Newton method
can be understood as a one-dimensional holomorphic dynamical system, and this was used in \cite{NewtonHSS} to give a good set of starting points that are guaranteed to find all roots. Other methods, such as the Ehrlich--Aberth method and the Weierstrass method, give rise to rational maps $\C^n \dasharrow \C^n$ for polynomials of degree $n$. In contrast to the Newton method, these higher-dimensional methods approximate all roots at the same time.

So far, not much is known about the global dynamics of higher-dimensional root-finding methods. In \cite{reinke_weierstrass_2020}, two phenomena were established for the Jacobi variant of the Weierstrass method: the existence of attracting cycles (and thereby non-general convergence), and the
existence of diverging orbits to infinity. This is a stark contrast to the Newton method: there $\infty \in \BP^1$ is a repelling fixed point.

In this paper we continue the investigation of diverging orbits. We show the following:
\begin{theorem}
  For every degree $d\geq 4$, there is a polynomial $\tilde{p}$ of degree $d$ with distinct roots such that the Jacobi variant of the Ehrlich--Aberth method for $\tilde{p}$ has diverging orbits in $\C^d$.
  \label{cor:eadeg}
\end{theorem}

\begin{theorem}
For every polynomial $\tilde{p}$ with distinct roots of degree $d \geq 3$, the Gauss--Seidel variant of the Weierstrass method for $\tilde{p}$ has diverging orbits in $\C^d$.
  \label{cor:wsdeg}
\end{theorem}

In Section~\ref{sec:hd} we provide common background of the discussed root-finding methods and a discussion for the various variants on how to use them.
We discuss the Jacobi variant of the Ehrlich--Aberth method in Section~\ref{sec:ea} and the Gauss--Seidel variant of the Weierstrass method in Section~\ref{sec:ws}.

\emph{Acknowledgements.} We gratefully acknowledge support by the Advanced Grant HOLOGRAM by the European Research Council.
The proof of both theorems were obtained with the help of computer algebra systems, in particular the \emph{Maxima}~\cite{Maxima} kernel of \emph{SageMath}~\cite{sagemath} and \emph{Singular}~\cite{DGPS}.

\section{Higher-dimensional root-finding methods}
\label{sec:hd}
In this section we define the maps under consideration. Let $p \in \C[z]$ be a monic polynomial of degree $n$. The root-finding methods
considered in this paper try to find all $n$ roots at the same time, by doing an iteration on $\C^n$. In fact, we
will consider the root-finding methods as rational maps $F \colon \C^n \dashrightarrow \C^n$, and try to extend them to sensible completions of $\C^n$.

We have an update function $g_p \colon \C \times \C^{n-1} \rightarrow \C$ where heuristically
$g_p(z_i;z_1,\dots, \hat{z_i},\dots z_n)$ tries to find a better approximation for $z_i$ to a root of $p$ under the assumption that all other coordinates
$z_1,\dots, \hat{z_i},\dots z_n$ are already close to roots of $p$. We will only consider update functions which are symmetric in the $\C^{n-1}$ component.

We consider two kinds of using the update function $g_p$ on an approximation vector $(z_1,\dots,z_n)$. One possibility is to apply $g_p$ simultaneously on every component, this is called  the Jacobi variant. Since we assume that $g_p$ is symmetric in the $\C^{n-1}$ component, the Jacobi variant is $S_n$ equivariant.

Another possibility is to apply the update function one at a time, and using the results already for the coming coordinates. This is called the Gauss--Seidel variant. See the algorithm listings. Our naming scheme comes from classical methods in numerical linear algebra.%

\begin{algorithm}
  \caption{Jacobi variant}
  \begin{algorithmic}
   \Function{Jacobi}{$z_1,\dots,z_n$}
   \For{$i \in 1,\dots,n$}
   \State $z'_i \leftarrow g_p(z_i;z_1,\dots, \hat{z_i},\dots z_n)$
   \EndFor
    \State return $(z'_1,\dots,z'_n)$
   \EndFunction
  \end{algorithmic}
  \label{algo:j}
\end{algorithm}

\begin{algorithm}
  \caption{Gauss--Seidel variant}
  \begin{algorithmic}
    \Function{GaussSeidel}{$z_1,\dots,z_n$}
    \For{$i \in 1,\dots,n$}
    \State $z_i \leftarrow g_p(z_i;z_1,\dots, \hat{z_i},\dots z_n)$
   \EndFor
    \State return $(z_1,\dots,z_n)$
   \EndFunction
  \end{algorithmic}
  \label{algo:gs}
\end{algorithm}

Both the Weierstrass method and Ehrlich--Aberth method use the rational function

\[q_i(z) = \frac{p(z)}{\prod_{j\not=i}(z-z_j)}\] as basis for their heuristic step. If $p(z) = \prod_j (z - \alpha_j)$ and all $z_j$ apart from $z_i$ are ``close'' to roots $\alpha_j$ of $p$, then
  $q_i$ is ``close'' to the linear function $(z - \alpha_i)$. The Weierstrass update $\operatorname{W}_p(z_i;z_1,\dots, \hat{z_i},\dots z_n)$ tries to solve $q_i$, pretending that it is a linear function, so
  \begin{equation}
    \operatorname{W}_p(z_i;z_1,\dots, \hat{z_i},\dots z_n) = z_i - q_i(z_i)\text{.}
    \label{def:wsstep}
  \end{equation}

  For the Ehrlich--Aberth function, we do a Newton step

  \begin{equation}
    \operatorname{EA}_p(z_i;z_1,\dots, \hat{z_i},\dots z_n) = z_i - \frac{q_i(z_i)}{{q'}_i(z_i)}\text{.}
    \label{def:eastep}
  \end{equation}
\begin{remark}
  We can also think of the Gauss--Seidel variant as the $n$-th iterate of the function
  $R_p \colon (z_1,\dots,z_n) \mapsto (z_2,\dots,z_n,g_p(z_1;z_2,\dots,z_n))$. This is a kind of cyclic shift,
  it takes the first coordinate, applies the step function, and moves it to the end. Since we assume that the step functions
  are symmetric in the additional coordinates, it is easy to see that $R^n_p$ is indeed the whole Gauss-Seidel step.
  \label{rem:cyclic}
  \end{remark}
\section{Ehrlich--Aberth}
\label{sec:ea}
\subsection{Möbius invariance}

\begin{lemma}[Field interpretation]
  Let $p(z) = \prod_i (z - \alpha_i)$ and $(z_1,\dots,z_n) \in \C^n$.
  Generically,
  \begin{equation}
    z_i' = \operatorname{EA}_p(z_i;z_1,\dots, \hat{z_i},\dots z_n) \Leftrightarrow \frac{1}{z_i - z_i'} = \sum_j \frac{1}{z_i - \alpha_j} - \sum_{j\not=i} \frac{1}{z_i - z_j}
    \label{eqn:field}
  \end{equation}
  \label{lem:field}
\end{lemma}
\begin{proof}
  Note that $q_i / q'_i$ is the inverse of the logarithmic derivative of $q_i$, so
  \[
    q_i(z_i) / q'_i(z_i) = \frac{1}{\sum_j \frac{1}{z_i - \alpha_j} - \sum_{j\not=i} \frac{1}{z_i - z_j}}\text{.}
  \]
  The claim now follows from elementary algebra.
\end{proof}
Without the second sum in the right hand side of \eqref{eqn:field}, we would obtain the update step for the classical Newton method. We can think of the Ehrlich--Aberth method as a correction of the Newton method, where we make add a term to avoid that coordinates converge to the same root. See Aberth's derivation~\cite{Aberth1973} for an electrostatic interpretation of Lemma~\ref{lem:field}.
\begin{lemma}[Möbius invariance]
  Let $M(z) = \frac{az+b}{cz+d}$ be a Möbius transformation, let $p(z) = \prod_i (z - \alpha_i)$,
  let $\tilde{p}(z) = \prod_i (z - M(\alpha_i))$,
  then $M(\operatorname{EA}_p(z_i;z_1,\dots, \hat{z_i},\dots,z_n)) = \operatorname{EA}_{\tilde{p}}(M(z_i);M(z_1),\dots, \hat{M(z_i)},\dots M(z_n))$
as rational functions in $z_i;z_1,\dots, \hat{z_i},\dots, z_n$.
  \label{lem:moebinv}
\end{lemma}
\begin{proof}
  By the previous lemma, it is enough to show that
  \[
  \frac{1}{z_i - z_i'} = \sum_j \frac{1}{z_i - \alpha_j} - \sum_{j\not=i} \frac{1}{z_i - z_j} \Leftrightarrow
  \frac{1}{M(z_i) - M(z_i')} = \sum_j \frac{1}{M(z_i) - M(\alpha_j)} - \sum_{j\not=i} \frac{1}{M(z_i) - M(z_j)}
\]
holds generically. This is a straightforward computation for a generating set of Möbius transformations such as translations, scalings and the involution $z \mapsto z^{-1}$.
\end{proof}
Motivated by the previous two lemmas, it makes sense to extend $\operatorname{EA}_p$ to a rational function $\operatorname{EA}_p \colon \PP^1 \times \left( \PP^1 \right)^{n-1} \rightarrow \PP^1$. One way to make this precise is to consider the ``universal Ehrlich--Aberth'' map: %
\begin{lemma}
  \label{lem:indet}
  The universal Ehrlich-Aberth map
  \begin{align*}
    \uEA \colon \PP^1 \times \left( \PP^1 \right)^{n-1} \times  \left( \PP^1 \right)^{n} &\dasharrow \PP^1 \\
    (z_1;z_2,\dots,z_n;\alpha_1,\dots,\alpha_n) &\mapsto EA_p (z_1;z_2,\dots,z_n) \text{ for } p(z) = \prod_i (z-\alpha_i)
  \end{align*}
  is Möbius equivariant. We have
  \begin{equation}
    \label{eqn:eainf}
    \uEA(\infty;w_2,\dots,w_n;\alpha_1,\dots,\alpha_n) = \sum^n_{i=1} \alpha_i - \sum^{n-1}_{j=1} w_j
  \end{equation}
  In particular, $\uEA(z_1;z_2,\dots,z_n;\alpha_1,\dots,\alpha_n) = z_1$ if $z_1$ is equal to exactly one of
  $z_2,\dots,z_n;\alpha_1,\dots,\alpha_n$, and the indeterminacy locus of $\uEA$ is contained in the
  union of intersections of two diagonals of the form $z_1 = z_j$ or $z_1 = \alpha_i$.
\end{lemma}
\begin{proof}
  Möbius equivariance is a reformulation of Lemma~\ref{lem:moebinv}.
  Equation~\eqref{eqn:eainf} then follows from equation~\eqref{eqn:field} by applying the Möbius transformation $z \mapsto z^{-1}$ for $z_1 = 0$.

  The right hand side of equation~\eqref{eqn:eainf} is well defined if at most one of the summands is equal to $\infty$.
  From this, and Möbius equivariance, the rest follows.
\end{proof}

One should note that we can now move some roots of the polynomials to $\infty$. In this case $\tilde{p}$ has lower degree then dimension $n$ of the dynamical space,
the formula \ref{def:eastep} still make sense in this case.

\subsection{Periodic Points of Permutation Type}
For a polynomial $p$ of degree $\leq n$, we denote the Jacobi variant using the Ehrlich--Aberth step associated to $p$ by $\operatorname{JEA}_p$.
So $\operatorname{JEA}_p \colon \left( \PP^1 \right)^{n} \dashrightarrow \left( \PP^1 \right)^{n}$.
By applying Lemma~\ref{lem:indet} componentwise, we see that the indeterminacy locus is contained in
the union of intersections of two diagonals of the form $z_i = z_j$ or $z_i = \alpha_j$.

  The easiest orbits to understand are the ones where one Jacobi step just
  permutes the coordinates. We will give an example with a nice geometric
  interpretation in the next section.
  However, we cannot hope to use them to find attracting cycles. Here is why:
  \begin{lemma}
    Let $p$ be a monic polynomial of degree $n$, let $z_1, \dots, z_n \in \PP^1$ be pairwise different.
    Assume none of the $z_i$ are roots of $p$ and that $\operatorname{JEA}_p(z_1, \dots , z_n) = \sigma_* (z_1, \dots , z_n)$ for some $\sigma \in \Sym(n)$.
    Then $\tr \TotalDiff{\sigma^{-1}_* \circ \operatorname{JEA}_p}{(z_1, \dots, z_n)} = -n$.
    In particular, the periodic orbit of $(z_{1}, \dots, z_n)$ is not attracting.
        \label{lem:trace}
  \end{lemma}
  \begin{proof}
    Since we assume that $z_1 \dots z_n$ are pairwise different and none of the $z_i$ are roots of $p$, the permutation $\sigma$ does not have a fixed point.
    By Möbius invariance, we can also assume that all $z_i$ are in $\C$. We claim that the matrix $\TotalDiff{\sigma^{-1}_* \circ \operatorname{JEA}_p}{(z_1, \dots, z_n)}$
    in the basis $\frac{\partial}{\partial z_1},\dots,\frac{\partial}{\partial z_n}$ has all diagonal entries equal to $-1$.
    For this it is enough to show the following: if $\operatorname{EA}_p(z_1;z_2,\dots,z_n) = z_2$, then $\frac{\partial \operatorname{EA}_p(z_1;z_2,\dots,z_n)}{\partial z_2} = -1$.
    This follows directly from Lemma~\ref{lem:field} by differentiating the field interpretation~\eqref{eqn:field} by $z_2$ and $z_1'$.

    In particular, the eigenvalues of $\TotalDiff{\sigma^{-1}_* \circ \operatorname{JEA}_p}{(z_1, \dots, z_n)}$
    can not all lie in $\D$. As the differential of the first return map for $(z_1, \dots, z_n)$ is a power of $\TotalDiff{\sigma^{-1}_* \circ \operatorname{JEA}_p}{(z_1, \dots, z_n)}$, the periodic orbit of $(z_{1}, \dots, z_n)$ is not attracting.
  \end{proof}
\subsection{Two-cycles via harmonic quadruples}
The idea in this section is to construct a number of periodic points of period 2 for generic degree 4 polynomials such that they have a nice description in terms of Möbius transformations.

We will use the following definition of the cross-ratio
\[\CR(a,b,c,d) = \frac{(a-c)(b-d)}{(a-d)(b-c)}\text{.}\]
The cross-ratio is well defined whenever $\left\{a,b,c,d\right\} \subset \PP^1$ has cardinality at least 3.
It is well known that the cross-ratio is invariant under Möbius transformations.
In our normalization $\CR(0,\infty,c,d) = c/d$.

For $a \not= b \in \PP^1$ let $M_{a,b} \in \PSL(2,\C)$ the unique Möbius involution which fixes $a$ and $b$.
It is easy to see that this is well-defined, for example $M_{0,\infty}(z) = -z$, so if $N \in \PSL(2,\C)$ moves
$a$ to $0$ and $b$ to $\infty$, then $M_{a,b} = N^{-1} \circ M_{0,\infty} \circ N$.
\begin{lemma}
  Let $a,  b, c, d \in \PP^1$ be pairwise different.
  The following are equivalent:
  \begin{enumerate}
    \item The cross-ratio $\CR(a,b,c,d)$ is equal to $-1$.
      \label{cr}
    \item $M_{a,b}$ interchanges $c$ and $d$.
      \label{interchange}
    \item $M_{c,d}$ interchanges $a$ and $b$.
  \end{enumerate}
We say in this case
that $a,b,c,d$ are harmonic.
\label{lem:harmonic}
\end{lemma}
\begin{proof}
  Since $\CR(a,b,c,d) = \CR(c,d,a,b)$, it is enough to show that (\ref{cr}) and (\ref{interchange}) are equivalent.
  By Möbius invariance, we can assume wlog.\ that $a = 0$ and $b = \infty$. In this case $M_{0,\infty}(z) = -z$,
  and $\CR(0,\infty,c,d) = c/d$, so the statement is obvious.
\end{proof}
\begin{remark}
  If $a,b,c,d$ are harmonic, they must lie on a real circle or a real line. If $a,b,c,d$ lie on a
real line, this definition extends the classical notion of harmonic quadruples. It is
easy to check that this is really a property of the unordered pair of unordered
pairs $\left\{ \left\{ a,b \right\},\left\{ c,d \right\} \right\}$.

Since cross-ratios are invariant under Möbius transformations, so is the notion of harmonic quadruples.
  \label{rem:harmonic}
\end{remark}
\begin{lemma}
  Let $a,  b, c, d \in \PP^1$ be pairwise different. Then there is a unique unordered pair $\left\{ e,f \right\}$ such that $a,b$ is harmonic to $e,f$ and $c,d$ is harmonic to $e,f$.  \label{lem:harmonicint}
\end{lemma}
\begin{proof}
  Using criterion (\ref{interchange}) of Lemma~\ref{lem:harmonic}, it is straightforward to verify that $e,f$ are the two fixed points of $M_{a,b} \circ M_{c,d}$.
\end{proof}
\begin{proposition}
  Let $p(z) = \prod^4_{i=1}(z-\alpha_i)$ be a polynomial with distinct roots. Let $(z_1,z_2,z_3,z_4)$ be such that
  \begin{itemize}
    \item
      $\left\{ z_1,z_2 \right\}$ is the unique
      unordered pair harmonic to $\left\{ \alpha_1, \alpha_2 \right\}$ and $\left\{ \alpha_3, \alpha_4 \right\}$,
     \item $\left\{ z_3, z_4 \right\}$
      is the unique unordered pair harmonic to $\left\{ \alpha_1, \alpha_3 \right\}$ and $\left\{ \alpha_2, \alpha_4 \right\}$.
  \end{itemize}
  Then under the Jacobi variant of the Ehrlich--Aberth method, $(z_1, z_2, z_3, z_4)$ is sent to $(z_2, z_1, z_4, z_3)$.
  \label{thm:geometricintr}
\end{proposition}
\begin{proof}
  By symmetry, it is enough to show that $\operatorname{EA}_p(z_1;z_2,z_3,z_4) = z_2$. By Möbius invariance, we can assume that $z_1 = \infty$, $z_2 = 0$, and
  in this case $\operatorname{EA}_p(\infty; 0;z_3,z_4) = (\sum^4_{i=1} \alpha_i) - z_3 - z_4$. Now $\left\{ \infty, 0 \right\}$ is harmonic to $\left\{ \alpha_1, \alpha_2 \right\}$ and harmonic to $\left\{ \alpha_3, \alpha_4 \right\}$, so $M_{0,\infty}(z) = -z$ both interchanges $\alpha_1$ and $\alpha_2$ as well as $\alpha_3$ and $\alpha_4$. Hence $M_{0,\infty}(z)$ fixes the unordered pair of unordered pairs $\left\{ \left\{ \alpha_1, \alpha_3 \right\}, \left\{ \alpha_2, \alpha_2 \right\}  \right\}$ set-wise, so  $M_{0,\infty}$ must also interchange $z_3$ and $z_4$. So $\alpha_1 = -\alpha_2, \alpha_3 = - \alpha_4, z_3 = -z_4$ and thus
  $\operatorname{EA}_p(\infty; 0,z_3,z_4) = 0$.
\end{proof}
\begin{remark}
  Using this proposition, we can produce $\binom{4}{2}\cdot2\cdot2 = 24$ periodic points for the Jacobi variant of the Ehrlich--Aberth method with mapping behavior
  $(z_1, z_2, z_3, z_4) \mapsto (z_2, z_1, z_4, z_3)$. These are far from the only ones, computer algebra shows that there are generically 72 such periodic points. We are not aware of a description of the 48 remaining ones that is as nice as given in this theorem.
\end{remark}
\subsection{Diverging Orbits for Jacobi Ehrlich--Aberth}
By Proposition~\ref{thm:geometricintr}, we have a nice class of periodic orbits for the Jacobi Ehrlich--Aberth method. Here is a particularly nice parametrization:
\begin{lemma}
  Let $\lambda \in \C \setminus \left\{ -1, 0, 1 \right\}$. For $p_\lambda(z) = (z^2-1)(z^2-\lambda^4)$, the points $(\lambda,-\lambda,\sqrt{-1}\lambda,-\sqrt{-1}\lambda)$
  and $(\lambda,-\lambda,0,\infty)$ are periodic points of period 2 with mapping scheme $(z_1, z_2, z_3, z_4) \mapsto (z_2, z_1, z_4, z_3)$.
  \label{lem:param}
\end{lemma}
\begin{proof}
  Straight forward computations using Lemma~\ref{lem:harmonic} show that
  $\left\{ \lambda, -\lambda \right\}$ is harmonic to $\left\{ 1,\lambda^2
  \right\}$ and harmonic to $\left\{ -1,-\lambda^2 \right\}$; that $\left\{
  \sqrt{-1}\lambda, -\sqrt{-1}\lambda \right\}$ is harmonic to $\left\{ -1,
  \lambda^2 \right\}$ and harmonic to $\left\{ 1, -\lambda^2 \right\}$; and
  that $\left\{ 0,\infty \right\}$ is harmonic to $\left\{ 1,-1 \right\}$ and harmonic to $\left\{ \lambda^2,-\lambda^2 \right\}$.
  The claim now follows from Proposition~\ref{thm:geometricintr}.
\end{proof}
We can explicitly compute the first return maps for these periodic points.
\begin{proposition}[Two-cycles via harmonic configurations]
  Let $p$ be a polynomial of degree $4$ such that its roots form a parallelogram $\alpha_1,\dots,\alpha_4$.
  Then the Jacobi Ehrlich--Aberth method of degree $4$ for $p$ has a two-cycle in $(\PP^1)^4 \setminus \C^4$. Moreover, there is
  a nonempty open subset $U \subset \C$ such that if $\CR(\alpha_1,\alpha_2,\alpha_3,\alpha_4) \in U$, then this cycle has a stable manifold intersecting $\C^4$, hence the Jacobi Ehrlich--Aberth method has orbits diverging off to infinity.
  \label{thm:escorbits}
\end{proposition}
\begin{proof}
  The statement of the proposition is invariant under affine Möbius
  transformations. We can move any parallelogram to
  $(1,\lambda^2,-1,-\lambda^2)$ by an affine Möbius transformation, so it is
  enough to consider $p_\lambda(z) = (z^2-1)(z^2-\lambda^4)$ as in
  Lemma~\ref{lem:param}. In this case, we already know that
  $(\lambda,-\lambda,0,\infty)$ is part of a two-cycle in $(\PP^1)^4 \setminus
  \C^4$, so it remains to show for some open set $V \subset \C$, this cycle has
  a stable manifold intersecting $\C^4$ for $\lambda \in V$. The set $U$ will
  be the image of $V$ under $\lambda \mapsto CR(1,\lambda^2,-1,-\lambda^2)$.

  Let $K \colon (\PP^1)^4 \mapsto (\PP^1)^4$ be the permutation $(z_1, z_2, z_3, z_4) \mapsto (z_2, z_1, z_4, z_3)$. Then $(\lambda,-\lambda,0,\infty)$ is a fixed point of $K \circ \operatorname{EA}_{p_{\lambda}}$. We take as basis for the tangent space at the point $\partial z_1, \partial z_2, \partial z_3, z^{-2}_4 \partial z_4$, where $\partial z_i$ is the extension from $\C^4$ to $(\PP^1)^4$ of the constant vector field in direction of the $i$-th coordinate.
  In this basis, the differential $D(K \circ \operatorname{EA}_{p_{\lambda}})$ at $(\lambda,-\lambda,0,\infty)$ is given by
  \begin{equation}
    \left(\begin{array}{rrrr}
        -1 & -\frac{2 \, {\left({\lambda}^{4} + 14 \, {\lambda}^{2} + 1\right)}}{{\left({\lambda} + 1\right)}^{2} {\left({\lambda} - 1\right)}^{2}} & -4 & -4 \, {\lambda}^{2} \\
        -\frac{2 \, {\left({\lambda}^{4} + 14 \, {\lambda}^{2} + 1\right)}}{{\left({\lambda} + 1\right)}^{2} {\left({\lambda} - 1\right)}^{2}} & -1 & -4 & -4 \, {\lambda}^{2} \\
        -1 & -1 & -1 & -2 \, {\lambda}^{4} + 2 \, {\lambda}^{2} - 2 \\
        -\frac{1}{{\lambda}^{2}} & -\frac{1}{{\lambda}^{2}} & -\frac{2 \, {\left({\lambda}^{4} - {\lambda}^{2} + 1\right)}}{{\lambda}^{4}} & -1
    \end{array}\right)\text{.}
    \label{eqn:jacobian}
  \end{equation}

We use Sage to compute an explicit eigenspace decomposition:
we give a eigenspace decomposition of the matrix~\eqref{eqn:jacobian}, in dependence of $\lambda$.
Let $D = {{\lambda}^{16} - 8 \, {\lambda}^{14} + 12 \, {\lambda}^{12}
                + 8 \, {\lambda}^{10} + 230 \, {\lambda}^{8} + 8 \, {\lambda}^{6} + 12 \, {\lambda}^{4}
              - 8 \, {\lambda}^{2} + 1}$.
  This matrix has eigenvalues
  \begin{equation}
    \frac{{\lambda}^{4} + 30 \, {\lambda}^{2} + 1}{{\lambda}^{4} - 2 \, {\lambda}^{2} + 1}, -\frac{{\lambda}^{8} - {\lambda}^{6} + 16 \, {\lambda}^{4} - {\lambda}^{2} + \sqrt{D} + 1}{{\lambda}^{6} - 2 \, {\lambda}^{4} + {\lambda}^{2}}, -\frac{{\lambda}^{8} - {\lambda}^{6} + 16 \, {\lambda}^{4} - {\lambda}^{2} - \sqrt{D} + 1}{{\lambda}^{6} - 2 \, {\lambda}^{4} + {\lambda}^{2}}, \frac{2 \, {\lambda}^{4} - 3 \, {\lambda}^{2} + 2}{{\lambda}^{2}}
  \end{equation}
  with right eigenvectors given by the columns of the following matrix:
  \begin{equation}
\left(\begin{array}{rrrr}
1 & 1 & 1 & 0 \\
-1 & 1 & 1 & 0 \\
0 & \frac{{\lambda}^{8} - 4 \, {\lambda}^{6} - 10 \, {\lambda}^{4} - 4 \, {\lambda}^{2} + \sqrt{D} + 1}{8 \, {\left({\lambda}^{6} - 2 \, {\lambda}^{4} + {\lambda}^{2}\right)}} & \frac{{\lambda}^{8} - 4 \, {\lambda}^{6} - 10 \, {\lambda}^{4} - 4 \, {\lambda}^{2} - \sqrt{D} + 1}{8 \, {\left({\lambda}^{6} - 2 \, {\lambda}^{4} + {\lambda}^{2}\right)}} & 1 \\
0 & \frac{{\lambda}^{8} - 4 \, {\lambda}^{6} - 10 \, {\lambda}^{4} - 4 \, {\lambda}^{2} + \sqrt{D} + 1}{8 \, {\left({\lambda}^{8} - 2 \, {\lambda}^{6} + {\lambda}^{4}\right)}} & \frac{{\lambda}^{8} - 4 \, {\lambda}^{6} - 10 \, {\lambda}^{4} - 4 \, {\lambda}^{2} - \sqrt{D} + 1}{8 \, {\left({\lambda}^{8} - 2 \, {\lambda}^{6} + {\lambda}^{4}\right)}} & -\frac{1}{{\lambda}^{2}}
\end{array}\right)
  \end{equation}
  For $\lambda = \pm\frac{1}{2} \, \sqrt{\pm\sqrt{-7} + 3}$, we find that the eigenvalues
  are $-63, 61.2529526\dots, -2.2529526\dots$ and $0$. By small (but nonzero) perturbation of $\lambda$, we thereby have
  an open set $V \subset \C$ where the first three eigenvalues lie in $\C \setminus \overline \D$, and the last eigenvalue is in $\D^*$.
  In this case, the point $(\lambda,-\lambda,0,\infty)$ is hyperbolic for the  map $K \circ \operatorname{JEA}$ with one attracting direction,
  so by the stable manifold theorem (see for example \cite[Ch.~2, Section~6]{PalisdeMelo}), we have a complex one-dimensional manifold that
  is tangential to $(0,0,1,-\frac{1}{{\lambda}^{2}})$, so in a neighborhood, points on the stable manifold give rise to diverging orbits in $\C^4$.
\end{proof}
\subsection{Extension to higher degrees}
We want to extend this result to higher degrees. We can use the following embedding:
\begin{lemma}
 Let $p$ be a polynomial of degree $d$, let $\alpha \in \C$ and let $\tilde{p}(z) =(z - \alpha)p(z)$.
We have commutative diagrams of rational maps
      \begin{tikzcd}[ampersand replacement=\&]
        (\BP^1)^d \ar[d,"\iota_\alpha"] \ar[r,dashed,"\operatorname{JEA}_p"] \& (\BP^1)^d \ar[d,"\iota_\alpha"] \\
        (\BP^1)^{d+1} \ar[r,dashed,"\operatorname{JEA}_{\tilde p}"] \& (\BP^1)^{d+1}
      \end{tikzcd}
      where $ \iota_\alpha\left( z_1,\dots,z_n \right) = (z_1,\dots,z_n,\alpha)$. Moreover, if $z_1,\dots,z_n \in \BP^1 \setminus \left\{ \alpha \right\}$,
      then $\operatorname{JEA}_p(z_1,\dots,z_n)$ is well defined if and only if $\operatorname{JEA}_{\tilde {p}}(z_1,\dots,z_n,\alpha)$ is.
  \label{ea:degemb}
\end{lemma}
\begin{proof}
  This follows directly from Lemma~\ref{lem:field}.
\end{proof}
\begin{proof}[Proof of Theorem \ref{cor:eadeg}]
  Let $k = d - 4$.
  Let $\lambda \in U$ as in Proposition~\ref{thm:escorbits}. Let $\alpha_1, \dots, \alpha_k \in \C \setminus \left\{ 0, \lambda, -\lambda \right\}$ and
  take $\tilde{p} = \prod_i (z-\alpha_i) p(z)$. We use a repeated application of Lemma~\ref{ea:degemb} to embed the dynamics of $\operatorname{JEA}_p$ into
  $\operatorname{JEA}_{\tilde{p}}$ via the map $(z_1,\dots,z_4) \mapsto (z_1,\dots,z_4,\alpha_1,\dots,\alpha_k)$. In a neighborhood of the periodic point $(0,\infty,\lambda,-\lambda)$, we see that $\operatorname{JEA}_{\tilde{p}}$ is well defined on the image of the stable manifold for $(0,\infty,\lambda,-\lambda)$, and this again gives rise to diverging orbits in $\C^d$. %
\end{proof}
\section{Weierstrass}
For a polynomial $p$ of degree $\leq n$, we denote the Gauss--Seidel variant using the Weierstrass step associated to $p$ by $\operatorname{GSW}_p$,
and the cyclic shift from Remark~\ref{rem:cyclic} for the Weierstrass map for cubic polynomials by $R_p$, so
\begin{eqnarray*}
  R_p(z_1,z_2,z_3) = \left(z_2,z_3,z_1 - \frac{p(z_1)}{(z_1-z_2)(z_1-z_3)}\right)\text{.}
\end{eqnarray*}
\label{sec:ws}
\subsection{Affine invariance}
\begin{lemma}[Affine invariance]
  Let $M(z) = {az+b}$ be an affine transformation let $p(z) = \prod_i (z - \alpha_i)$,
  let $\tilde{p}(z) = \prod_i (z - M(\alpha_i))$,
  then $M(\operatorname{W}_p(z_i;z_1,\dots, \hat{z_i},\dots z_n)) = \operatorname{W}_{\tilde{p}}(M(z_i);M(z_1),\dots, \hat{M(z_i)},\dots M(z_n))$
  as rational functions in $(z_i;z_1,\dots, \hat{z_i},\dots z_n)$
  \label{lem:affinv}
\end{lemma}
\begin{proof}
  Straight forward computation.
\end{proof}
\subsection{Special case $z^3$}
For $p=z^3$, the Gauss--Seidel variant of the Weierstrass map is homogeneous, in the sense that $\operatorname{GSW}_{z^3}(\lambda z_1, \lambda z_2, \lambda z_3) = \lambda \operatorname{GSW}_{z^3}(z_1, z_2, z_3)$. This is also true for the rotational map $F_{z^3}$ from Remark~\ref{rem:cyclic}. So the map $F_{z^3}$ fibers as rational map over $\phi \colon \BP^2 \dasharrow \BP^2$, where
\begin{eqnarray*}
  \phi( [z_1 : z_2 : z_3 ]) = [z_2 : z_3 : z_1 - \frac{z_1^3}{(z_1-z_2)(z_1-z_3)}] = \big[&(z_1-z_2)(z_1-z_3)z_2 \\
    :&(z_1-z_2)(z_1-z_3)z_3 \\
  :&(z_1-z_2)(z_1-z_3)z_1 - z_1^3\big] \text{.}
\end{eqnarray*}
Since $F^3_p = \operatorname{GSW}_p$, the Gauss-Seidel variant of the Weierstrass map for $z^3$ fibers over $\phi^3$.
\begin{lemma}
  The map $\phi$ has a periodic point $[z_1 : z_2 : z_3]$ of period $3$ such that no coordinate
  vanishes, and $\operatorname{GSW}_{z^3}(z_1,z_2,z_3) = \mu (z_1, z_2, z_3)$ for some $\mu \in \C$ with $|\mu| > 1$.
  Moreover, we can choose a periodic point such that the cycle is repelling for $\phi \colon \BP^2 \dasharrow \BP^2$.
  \label{lem:z3esc}
\end{lemma}
\begin{proof}
  We establish this using explicit computations for the map $\operatorname{GSW}_{z^3}$. Let us use variables $z_1,\dots,z_6$. We impose two relations: $(z_4,z_5,z_6)$ should be the
  image of $(z_1,z_2,z_3)$ under the map $\operatorname{GSW}_{z^3}$, and also $(z_4, z_5, z_6)$ should be a scalar multiple of $(z_1,z_2,z_3)$.
  To avoid points of indeterminacy, it is more convenient to impose that $(z_1,z_2,z_3)$ is a scalar multiple of $(z_4,z_5,z_6)$.
  For this we introduce another variable $\lambda$, which will be the inverse of $\mu$.
  Away from points of indeterminacy, the relation $\operatorname{GSW}_{z^3}(z_1,z_2,z_3) = (z_4, z_5, z_6)$ is equivalent to the set of equations
  \begin{eqnarray*}
    (z_1 - z_2)(z_1 - z_3)(z_1 - z_4) - z_1^3  = 0\\
      (z_2 - z_3)(z_2 - z_4)(z_2 - z_5) - z_2^3  = 0\\
      (z_3 - z_4)(z_3 - z_5)(z_3 - z_6) - z_3^3  = 0 \text{.}\\
  \end{eqnarray*}
  The relation $\lambda(z_4,z_5,z_6) = (z_1,z_2,z_3)$ is equivalent to
  \begin{eqnarray*}
      -\lambda z_4 + z_1 = 0 \\
      -\lambda z_5 + z_2  = 0\\
      -\lambda z_6 + z_3 = 0 \text{.}
  \end{eqnarray*}
  We know that the space of possible solutions is homogeneous, so we introduce an extra normalization $z_6 = 1$. Now
  we have a polynomial system of equations that is easily solvable with \emph{Singular} or \emph{Maxima}.
  \emph{Singular} shows this has 18 solutions, and $\lambda$ has to be a root of
  $(\lambda^5 - 2\lambda^4 - 3\lambda^3 + 8\lambda^2 - 4\lambda + 1)(\lambda^3 - 5\lambda^2 + 4\lambda - 1)(\lambda^2 - 3\lambda + 1)$.
  For the factor $(\lambda^2 - 3\lambda +1)$, some coordinates of the cycles are $0$. For the other two factors, we can find solutions
  with $\lambda \in \D^*$, such that all $z_1,\dots,z_6$ are non-zero, and compute the differential of $\phi$ so see that
  there are also solutions that are repelling for $\phi$. %
\end{proof}
\subsection{Diverging Orbits for Gauss--Seidel Weierstrass variant}
Let $p$ be a monic cubic polynomial.
We extend $R_p$ to a map on $\PP^3$. Let $p^*(z,w)$ be the homogenization of $p$, i.e. $p^*(z,w) = p(z/w)w^3$.
Then $R_p$ has the following extension:
\begin{eqnarray*}
  R_p([z_1\colon z_2 \colon z_3 \colon w]) = \big[&(z_1-z_2)(z_1-z_3)z_2 \\
    :&(z_1-z_2)(z_1-z_3)z_3 \\
  :&(z_1-z_2)(z_1-z_3)z_1 - p^*(z_1,w) \\
  :& (z_1-z_2)(z_1-z_3)w\big]
  \label{eqn:R_p_homo}
\end{eqnarray*}
As Gauss--Seidel variant of the Weierstrass map is the third iterate of $R_p$, we focus on
$R_p$ for now.
\begin{lemma} %
  The restriction of $R_p$ to the plane $w = 0$ is independent of $p$ and is given by $\phi$. Moreover,
  if $[z_1:z_2:z_3]$ is periodic for $\phi$ of period $n$, with $F^n_{z^3}(z_1,z_2,z_3) = \mu (z_1,z_2,z_3)$,
  then the eigenvector of $D(F^n_p)$ at $[z_1:z_2:z_3:0]$ transverse to the plane $w = 0$ has eigenvalue $\mu^{-1}$.
  \label{lem:triangluar}
\end{lemma}
\begin{proof}
  On the plane $w = 0$, the Jacobian matrix associated to $R_p$ by equation~\eqref{eqn:R_p_homo} for has a $3 + 1$ triangular block decomposition,
  where the first block is the Jacobian matrix for $\phi$, and the last diagonal entry is independent of $p$.
  From this we see that the eigenvalues of $D(R^n_p)$ for a periodic point $[z_1:z_2:z_3:0]$ are independent of $p$. For $p = z^3$, it is clear that
  the transverse eigenvector has eigenvalue $\mu^{-1}$.
\end{proof}
\begin{proposition}
  For cubic polynomials $p$ with distinct roots, the Gauss--Seidel Weierstrass variant has diverging orbits that go to infinity in every component.
  \label{thm:gsws}
\end{proposition}
\begin{proof}
  In the language of Lemma~\ref{lem:triangluar}, the periodic point constructed in Lemma~\ref{lem:z3esc}
  is a periodic point for $R_p$ on the plane $w=0$ that has an attracting direction transverse to the plane, is repelling on $w = 0$,
  and has no other vanishing coordinate.
  Now we can again apply the stable manifold theorem to see that the points in the stable manifold of the periodic point are the diverging orbits we want.
\end{proof}
\subsection{Extension to higher degrees}

\begin{lemma}
  Let $p$ be a polynomial of degree $d$, let $\alpha \in \C$ and let $\tilde{p}(z) =(z - \alpha)p(z)$. For the Gauss--Seidel variant of the Weierstrass method,
we have commutative diagrams of rational maps
      \begin{tikzcd}[ampersand replacement=\&]
        \C^d \ar[d,"\iota_\alpha"] \ar[r,dashed,"\operatorname{GSW}_p"] \& \C^d \ar[d,"\iota_\alpha"] \\
        \C^{d+1} \ar[r,dashed,"\operatorname{GSW}_{\tilde p}"] \& \C^{d+1}
      \end{tikzcd}
      where $ \iota_\alpha\left( z_1,\dots,z_n \right) = (z_1,\dots,z_n,\alpha)$. Moreover, if $z_1,\dots,z_n \in \C \setminus \left\{ \alpha \right\}$,
      and $\operatorname{GSW}_p(z_1,\dots,z_n)$ is well defined and has no coordinate equal to $\alpha$, then $\operatorname{GSW}_{\tilde {p}}(z_1,\dots,z_n,\alpha)$ is well defined and
      equal to $\iota_\alpha\left( z_1,\dots,z_n \right)$.
  \label{lem:degembws}
\end{lemma}
\begin{proof}
  This follows directly from the defining equation~\ref{def:wsstep}, by stepping through the algorithm~\ref{algo:gs}.
  A key point here is that for the algorithm applied for $\operatorname{GSW}_{\tilde{p}}(z_1,\dots,z_n,\alpha)$, the contribution of $\alpha$ cancels under the assumption $z_i \not= \alpha$, the final step is $W_{\tilde{p}}(\alpha; z'_1,\dots,z'_n)$, there the second part of the assumption is used.
\end{proof}
\begin{proof}[Proof of Theorem~\ref{cor:wsdeg}]
  Let $k = d - 3$ and choose a factorization $\tilde{p}(z) = p(z) \prod_i(z - \alpha_i)$, where $\alpha_1,\dots,\alpha_k \in \C$ are pairwise different.
  By Proposition~\ref{thm:gsws}, there are $z^j_i \in C, 1\leq i \leq 3, j \in \N$, such that $(z^j_1, z^j_2, z^j_3)$ is an infinite forward orbit under the Gauss--Seidel variant for $p$, with $z^j_i \rightarrow \infty$ for all $1 \leq i \leq 3$. By taking a tail of the orbit, we can assume that none of the coordinates $z^j_i$ are equal to roots of $\tilde{p}(z)$. Now we can iteratively apply Lemma~\ref{lem:degembws} to see that $\left( z^j_1, z^j_2, z^j_3, \alpha_1,\dots,\alpha_k \right)$ is a well defined forward orbit for the Gauss--Seidel variant of the Weierstrass method for $\tilde{p}$,
  and the first 3 components go off to infinity.
  \end{proof}
\section{Outlook}
In this paper, we extended the results about diverging orbits for the Weierstrass method by changing
the step function to the Ehrlich--Aberth method or by changing the iteration procedure to the Gauss--Seidel variant.
A natural next step is to consider the Gauss--Seidel variant of the Ehrlich--Aberth method. Based on the Möbius invariance of
the Ehrlich--Aberth method, this is equivalent to finding periodic cycles with at least one attracting direction.
However, we are not aware of explicit descriptions of a family of
periodic points as in the Jacobi case.

Also, the existence of attracting cycles is still open for the various methods discussed here apart from the Jacobi variant of the Weierstrass method.
An intermediate question would be the existence of periodic points with multiple attracting eigenvalues.
\bibliographystyle{alpha}