\section{Introduction}

The fading channel is a widely adopted time-varying model for real-world wireless communications. In this model, the channel gain changes over time satisfying a certain distribution, called the channel distribution information (CDI). In a block fading channel model, the channel gain is assumed to be constant over a fixed time interval $T_c$, known as the coherent time, and change to a new independent value afterwards. In many of today's communication systems, channel estimation is performed in the first place to obtain the instantaneous channel state information (CSI), and then data transmission follows. However, in many scenarios, the coherence time is very short (e.g., only a few symbol intervals). In this case, channel estimation may significantly reduces the overall data rate. Besides, the estimation precision is quite limited. Consequently, communication without instantaneous CSI (or noncoherent communication) is preferable.

Polar codes are the first family of codes that provably achieves the capacity of any binary-input symmetric memoryless channels with low encoding and decoding complexity \cite{arikan2009channel}. Later, polar codes are generalized to asymmetric channels while still capacity-achieving \cite{honda2013asymmetric}. There have been studies on polar coding for fading channels under various CSI assumptions. In \cite{boutros2013fading}, polar coding for quasi-static fading channels with two states was studied. Polar coding for block fading channels with full CSI and i.i.d. fading channels with CDI was considered in \cite{Santos2013rayleigh}. For block fading binary symmetric and
additive exponential noise channels with CSI at the receiver (CSI-R), a hierarchical polar coding scheme was proposed in \cite{si2014fading}, which achieves capacity, but only works for block fading channels with finite states. A simple method for construction of polar codes for Rayleigh fading channel was presented in \cite{trifonov2015rayleigh}. Polar codes and polar
lattices for i.i.d. fading channels with CSI-R were constructed in \cite{liu2016fading}, which achieve the ergodic capacity through single-stage polarization. All of the aforementioned polar coding schemes for block fading channels require the coherent time to be very large. As far as we know, polar coding for block fading channels with arbitrary finite coherent time has not been investigated in literature yet.

By viewing transmitted symbols in a coherent block as a supersymbol, we can design polar codes for block fading channels by using techniques in coded modulation, i.e., multilevel coding (MLC) \cite{Imai1977mlc,Wachsmann1999mlc} and bit-interleaved coded modulation (BICM) \cite{Caire1998bicm,Guilleni2008bicm}. Polar coded modulation, with both the MLC approach and the BICM approach, has been investigated, e.g., \cite{Seidl2013pcm,Mahdavifar2016bicm}. It is shown that the MLC-based polar coded modulation scheme achieves the code modulation capacity, while the BICM-based scheme suffers a certain rate loss.

In this paper, we aim to design capacity-achieving polar codes for block fading channels, especially for those with only CDI. A block fading channel with coherent time $T_c$ can be decomposed into $T_c$ parallel sub-channels. In the case when only CDI is known by the communicators, the sub-channels are correlated, and we use an MLC-based approach to design polar codes since such an approach preserves channel capacity. When instantaneous CSI is available at the receiver side or at both sides, the sub-channels become independent from each other. Thus, we can use an interleaver to randomize symbols from different sub-channels, and perform joint encoding and decoding so as to enhance the finite length performance. We refer to this scheme as the BICM-based approach in this paper. We show that in all the above mentioned cases, our proposed schemes achieve the ergodic capacity of any binary-input memoryless block fading channel.

The rest of this paper is organized as follows. In Section \ref{S:problemstate} we introduce the block fading channel model and describe the main idea of our schemes. Section \ref{S:PolarPri} provides some related knowledge on polar codes. Details of our proposed schemes are presented in Section \ref{S:scheme}. Section \ref{S:Numerical} shows some numerical results on achievable rates of block Rayleigh fading channels and the synthesized sub-channels in our scheme. Section \ref{S:discuss} concludes this paper with some discussions.

\section{Problem Statement}
\label{S:problemstate}

We consider a block fading channel model with coherent interval $T_c$. At time interval $i$ ($i=1,2,...$), the channel is modeled as
\begin{equation}
\mathbf{y}^i=h^i\mathbf{x}^i+\mathbf{w}^i,
\end{equation}
where $h^i\in \mathbb{R}$ is the channel gain at time interval $i$, $\mathbf{x}^i=[x^i_1,x^i_2,...,x^i_{T_c}]^{\mathrm{T}}\in\{-1,1\}^{T_c}$ is the binary input signal after BPSK modulation, $\mathbf{y}^i=[y^i_1,y^i_2,...,y^i_{T_c}]^{\mathrm{T}}\in \mathbb{R}^{T_c}$ is the channel output, and $\mathbf{w}^i=[w^i_1,w^i_2,...,w^i_{T_c}]^{\mathrm{T}}\in \mathbb{R}^{T_c}$ is the white Gaussian noise, with $w^i_j\sim \mathcal{N}(0,\sigma^2)$ for every $j\in[T_c]$. Let $X_j$ and $Y_j$ ($j\in[T_c]$) respectively be the random variables standing for the $j$th input and output symbols of a coherent block, and $H$ be the random variable for the channel gain.

First, let us discuss the case when both the transmitter and the receiver only have the CDI of the channel. Consider a series of transmissions over $N$ fading blocks. In this paper, we call the $N$ consecutive coded blocks a \textit{frame}. Denote $\mathbf{X}=[\mathbf{x}^1,...,\mathbf{x}^N]$ and  $\mathbf{Y}=[\mathbf{y}^1,...,\mathbf{y}^N]$, and let $\mathbf{z}_j$ ($j\in[T_c]$) denote the $j$th row vector of $\mathbf{X}$. Then the mutual information of a transmission frame can be expanded as
\begin{align}
I(\mathbf{X};\mathbf{Y}
&=\sum_{j=1}^{T_c}I(\mathbf{z}_j;\mathbf{Y}|\mathbf{z}_{1:j-1}),\label{expansion}
\end{align}
where $\mathbf{z}_{1:j-1}$ is short for $\{\mathbf{z}_1,...,\mathbf{z}_{j-1}\}$. Similar abbreviations will be used throughout this paper. Note that
\begin{align}
\lim\limits_{N\rightarrow\infty}\frac{1}{N}I(\mathbf{z}_j;\mathbf{Y}|\mathbf{z}_{1:j-1})&=I(X_j;Y_{1:T_c}|X_{1:j-1})\nonumber\\
&=I(X_j;Y_{1:T_c},X_{1:j-1}),\nonumber
\end{align}
which is the mutual information of a binary-input channel
\begin{equation}
W^{(j)}(\mathbf{y},x_{1:j-1}|x_j):\{-1,1\} \rightarrow \mathbb{R}^{T_c}\times\{-1,1\}^{j-1}, \label{subchannel}
\end{equation}
with the $j$th input symbol of a block being the input, and the whole output $\mathbf{y}$ together with the previous $j-1$ input symbols of the block being the output.

Based on this expansion, we can use an MLC-based approach to design polar codes for block fading channels with only CDI. The encoding of a frame consists of $T_c$ component polar codes, designed for each of the $T_c$ sub-channels respectively. When an encoded frame is generated, the sender transmits it block by block. Having received a signal frame, the receiver uses a multistage decoder to decode the component polar codes one by one. At stage $j$ ($j\in[T_c]$), it decodes the $j$th sub-channel based on the received frame together with the estimates of previous stages. If the component polar codes are capacity-achieving, the ergodic CDI capacity of the binary-input block fading channel is also achievable with this scheme.

For the case when only the receiver knows or both the transmitter and the receiver know the instantaneous CSI, the conditional mutual information of a transmission frame can be expanded as
\begin{align}
I(\mathbf{X};\mathbf{Y}|\mathbf{h}
&=\sum_{j=1}^{T_c}I(\mathbf{z}_j;\mathbf{y}_j|\mathbf{h}),\label{CSIR-2}
\end{align}
where $\mathbf{h}=[h^1,...,h^N]^{\mathrm{T}}$ is the channel state vector for the $N$ blocks, and $\mathbf{y}_j$ is the $j$th row of $\mathbf{Y}$. Equation (\ref{CSIR-2}) can be proved by the fact that when $\mathbf{h}$ is known, $\mathbf{y}_j$ is only related to $\mathbf{z}_j$ and independent from $\mathbf{z}_i$ ($i\in[T_c]\setminus j$). Details of the proof are omitted in this paper due to space limitation.

Consider $\mathbf{z}_j$ and $\mathbf{y}_j$ respectively as the input and output of a synthesized channel $W^{(j)}$. Obviously $W^{(j)}$ is an i.i.d. fading channel with the same CDI as the block fading channel, the ergodic capacity of which under the CSI-R assumption is $I(X_j,Y_j|H)$, i.e.,  $\lim\limits_{N\rightarrow\infty}\frac{1}{N}I(\mathbf{z}_j;\mathbf{y}_j|\mathbf{h})=I(X_j,Y_j|H)$. Since $W^{(1)},...,W^{(T_c)}$ are equivalent to each other in this case, we can see that the block fading channel is equivalent to $T_c$ i.i.d. fading channels with identical instantaneous CSI. Thus, simply designing a CSI-R polar coding scheme \cite{liu2016fading}, which treats the channel gain as another channel output, for each of the sub-channels independently is sufficient to achieve capacity. Nevertheless, the transmitter can actually use an interleaver to merge the parallel channels into a single one so as to improve the finite length performance. We will refer to this approach as the BICM-based scheme. Suppose the transmitter randomly permutes the positions of symbols in a frame with an interleaver. By deinterleaving the received frame, the signals can be seen as being transmitted through an i.i.d. fading channel with the same CDI as the block fading channel \cite{tse2005fundamentals}, provided that the interleaver is perfect and $N$ is sufficiently large. Since the CDI of the merged i.i.d. fading channel is the same as the original one, the ergodic CSI-R capacity is also the same. Thus, a polar code designed for the merged channel can achieve the ergodic capacity of the block fading channel.

\section{Preliminaries on Polar Codes}
\label{S:PolarPri}

First, we fix some notations that will be used in the sequel. $N=2^n$ with $n$ being an arbitrary integer. $\mathbf{G}_N=\mathbf{B}_N \textbf{F}^{\otimes n}$ is the generator matrix of polar codes, with $\mathbf{B}_N$ being the bit-reversal matrix and $\textbf{F}=
\begin{bmatrix}
1 & 0 \\
1 & 1
\end{bmatrix}$. The Bhattacharyya parameter $Z(X|Y)$ of a random variable pair $(X,Y)$ is defined as
\begin{equation}
  Z(X|Y)=2\sum_{y\in \mathcal{Y}} P_Y(y)\sqrt{P_{X|Y}(0|y)P_{X|Y}(1|y)},
\end{equation}
with $X$ being binary and $Y$ being defined on an arbitrary discrete alphabet.

Let $X^{1:N}$ be $N$ independent copies of a binary random variable $X$, and $U^{1:N}=X^{1:N}\mathbf{G}_N$. It is shown that as $N$ goes to infinity, $U^{1:N}$ polarizes in the sense that $U^i$ ($i\in [N]$) becomes either almost independent of $U^{1:i-1}$ and uniformly distributed, or almost determined by $U^{1:i-1}$ \cite{arikan2010source}. Based on this phenomenon, for $\delta_N=2^{-N^\beta}$ with $\beta \in (0,1/2)$, we define the \textit{high entropy set} as
\begin{align}
\mathcal{H}^{(N)}_X&=\{i\in [N]:Z(U^i|U^{1:i-1})\geq 1-\delta_N\},\label{HX
\end{align}
which satisfies
\begin{equation}
\begin{aligned}
\lim_{N\rightarrow \infty}\frac{1}{N}|\mathcal{H}^{(N)}_X|&=H(X)
\end{aligned}
\end{equation}

Let $(X,Y)\sim p_{X,Y}$ be a random variable pair with $X$ being binary and $Y$ being defined on an arbitrary countable set. Consider $X$ as the source to be compressed and $Y$ as side information of $X$. Let $U^{1:N}=X^{1:N}\mathbf{G}_N$. Similar to the single source case, conditioned on $Y^{1:N}$, $U^{1:N}$ polarizes as $N$ goes to infinity. For $\delta_N=2^{-N^\beta}$ with $\beta \in (0,1/2)$, define the \textit{reliable set} as
\begin{align}
\mathcal{L}^{(N)}_{X|Y}&=\{i\in [N]:Z(U^i|Y^{1:N},U^{1:i-1})\leq \delta_N\},\label{LXY}
\end{align}
which satisfy
\begin{equation}
\begin{aligned}
\lim_{N\rightarrow \infty}\frac{1}{N}|\mathcal{L}^{(N)}_{X|Y}|&=1-H(X|Y).
\end{aligned}
\end{equation}

Consider $N$ independent uses of a binary-input discrete memoryless channel $W(Y|X)$. Let $U^{1:N}=X^{1:N}\mathbf{G}_N$, and define $\mathcal{H}^{(N)}_X$ as in (\ref{HX}) and $\mathcal{L}^{(N)}_{X|Y}$ as in (\ref{LXY}). To construct a polar code for $W$, we define \cite{honda2013asymmetric}
\begin{align}
\mathcal{I}&\triangleq \mathcal{H}_X^{(N)}\cap \mathcal{L}_{X|Y}^{(N)}, \label{PC-I}\\
\mathcal{F}^r&\triangleq \mathcal{H}_X^{(N)}\cap (\mathcal{L}_{X|Y}^{(N)})^C, \label{PC-Fr}\\
\mathcal{F}^d&\triangleq (\mathcal{H}_X^{(N)})^C,\label{PC-Fd}
\end{align}
where $\mathcal{A}^C$ denotes the complement set of $\mathcal{A}$. The encoding is done by assigning $\{u^i\}_{i\in \mathcal{I}}$ with information bits and $\{u^i\}_{i\in \mathcal{F}^r}$ with uniformly distributed frozen bits (shared between the sender and the receiver), calculating $\{u^i\}_{i\in \mathcal{F}^d}$ with
\begin{equation}
\label{DetBit}
u^i= \arg\max_{u=\{0,1\}} P_{U^i|U^{1:i-1}}(u|u^{1:i-1}),
\end{equation}
and finally computing $x^{1:N}=u^{1:N}\mathbf{G}_N$ since $\mathbf{G}_N=\mathbf{G}_N^{-1}$.

Upon receiving $y^{1:N}$, the receiver computes an estimate $\bar{u}^{1:N}$ of $u^{1:N}$ as
\begin{equation}
\bar{u}^{i}=
\begin{cases}
u^i,\text{~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~if } i\in \mathcal{F}^{r}\\
\arg\max_{u\in\{0,1\}}P_{U^{i}|U^{1:i-1}}(u|u^{1:i-1}),\text{~~~~if } i\in \mathcal{F}^{d}\\
\arg\max_{u\in\{0,1\}}P_{U^{i}|Y^{1:N}U^{1:i-1}}(u|y^{1:N},u^{1:i-1}),\\
\text{~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~if } i\in \mathcal{I}
\end{cases}.
\end{equation}

The rate of such a scheme, $R=|\mathcal{I}|/N$, satisfies
\begin{equation}
\lim_{N\rightarrow \infty}R=I(X;Y).
\end{equation}
The block error probability of such a scheme can be upper bounded by
\begin{equation}
\label{SC-EP}
P_e \leq\sum_{i\in\mathcal{I}}{Z(U^i|Y^{1:N},U^{1:i-1})}=O(2^{-N^\beta}).
\end{equation}

\section{Proposed Polar Coding Schemes}
\label{S:scheme}
\subsection{Only CDI Available}
The joint transition probability density function (PDF) of a coherent block without instantaneous CSI is given by
\begin{equation}
\label{transp}
p(\mathbf{y}|\mathbf{x})=\int_{0}^{\infty} \big{(}\prod_{j=1}^{T_c}p(y_j|h,x_j)\big{)}f(h) \mathrm{d}h,
\end{equation}
where $\mathbf{x}=(x_1,...,x_{T_c})$, $\mathbf{y}=(y_1,...,y_{T_c})$, $p(y_j|h,x_j)$ is the PDF for a given channel gain $h$, and $f(h)$ is the distribution of $h$. The transition PDF of the $j$th ($j\in[T_c]$) sub-channel is
\begin{equation}
p(\mathbf{y},x_{1:j-1}|x_j)=\sum_{x_{j+1:T_c}}p(x_{1:j-1},x_{j+1:T_c})p(\mathbf{y}|\mathbf{x}).\label{subchpdf}
\end{equation}

One may easily verify that the sub-channels are symmetric if the original fading channel is symmetric. In the symmetric case, uniform input distribution achieves the capacity, and the deterministic set defined in (\ref{PC-Fd}) is null. If the original channel is asymmetric, then the capacity-achieving input distribution may not be uniform distribution, and the deterministic set will not be empty.

For the $j$th ($j\in[T_c]$) sub-channel $W^{(j)}(\mathbf{y},x_{1:j-1}|x_j)$, let $\mathbf{u}_j=\mathbf{z}_j\mathbf{G}_N$. Define the high entropy set $\mathcal{H}^{(N)}_{X_j}$ in the same way as (\ref{HX}), and the reliable set $\mathcal{L}^{(N)}_{X_j|Y_{1:T_c}X_{1:j-1}}$ by
\begin{equation}
\begin{aligned}
\mathcal{L}^{(N)}_{X_j|Y_{1:T_c}X_{1:j-1}}&\triangleq\{i\in [N]:Z(U_j^i|Y_{1:T_c}^{1:N},U_{1:j-1}^{1:N},\\
&~~~~~~~~~~~~~~~~U_j^{1:i-1})\leq \delta_N\},
\end{aligned}
\end{equation}
with $U_{1:0}^{1:N}=\emptyset$. The multilevel encoding procedure goes as follows.
\begin{itemize}
	\item For the $j$th ($j\in[T_c]$) sub-channel, define the information bit set $\mathcal{I}_j$, frozen bit set $\mathcal{F}_j^r$ and deterministic bit set $\mathcal{F}_j^d$ according to (\ref{PC-I}), (\ref{PC-Fr}) and (\ref{PC-Fd}) respectively.
	\item Insert information bits to $\{u_j^i\}_{i\in \mathcal{I}_j}$ and frozen bits to $\{u_j^i\}_{i\in \mathcal{F}_j^r}$, and compute $\{u_j^i\}_{i\in \mathcal{F}_j^d}$ according to (\ref{DetBit}).
	\item Compute $\mathbf{z}_j=\mathbf{u}_j\mathbf{G}_N$ for each $j\in[T_c]$ and generate the final coded frame by $\mathbf{X}=\big{(}\mathbf{z}_1;...;\mathbf{z}_{T_c}\big{)}$.
\end{itemize}

The sender transmits $\mathbf{X}$ column by column. Having received $\mathbf{Y}$, the receiver performs multistage decoding.
In the $j$th ($1\leq j\leq T_c$) stage, the decoder decodes $\mathbf{u}_j$ with the aid of the estimates in previous stages:
	\begin{equation}
	\bar{u}_j^{i}=
	\begin{cases}
	u_j^{i},  \text{~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~if } i\in \mathcal{F}_j^r \\
	\arg\max_{u\in\{0,1\}}P_{U_j^{i}|U_j^{1:i-1}}(u|u_j^{1:i-1}),\text{~~if } i\in \mathcal{F}_j^d\\
	\arg\max_{u\in\{0,1\}}P_{U_j^{i}|Y^{1:N}_{1:T_c}U_{1:j-1}^{1:N}U_j^{1:i-1}}(u|y^{1:N}_{1:T_c},\\
	~~\bar{u}_{1:j-1}^{1:N},u_j^{1:i-1}),\text{~~~~~~~~~~~~~~~~~~~~~~~~~if } i\in \mathcal{I}_j
	\end{cases}
	\end{equation}
where $\bar{u}_{1:0}^{1:N}=\emptyset$.

The block error probability of the $j$th component polar code provided that the previous component codes are correctly decoded can be upper bounded by
\begin{equation}
P_e^{(j)} \leq\sum_{i\in\mathcal{I}_j}{Z(U_j^i|Y_{1:T_c}^{1:N},U_{1:j-1}^{1:N},U_j^{1:i-1})}=O(2^{-N^\beta})
\end{equation}
according to the definition of the information bit set. Thus, the overall frame error probability can be upper bounded by
\begin{equation}
\label{FER}
P_e \leq\sum_{j=1}^{T_c}P_e^{(j)}=T_cO(2^{-N^\beta}).
\end{equation}

The asymptotic rate of the $j$th ($j\in[T_c]$) component polar code is
\begin{equation}
\lim\limits_{N\rightarrow\infty}R_j=\lim\limits_{N\rightarrow\infty}\frac{|\mathcal{I}_j|}{N}=I(X_j;Y_{1:T_c}|X_{1:j-1}).
\end{equation}
Thus, the asymptotic rate of the scheme is
\begin{align}
\lim\limits_{N\rightarrow\infty}R&=\lim\limits_{N\rightarrow\infty}\frac{1}{T_c}\sum_{j=1}^{T_c}R_
=\frac{1}{T_c}I(X_{1:T_c};Y_{1:T_c}), \label{AR}
\end{align}
which equals the ergodic capacity of the block fading channel.

From (\ref{FER}) and (\ref{AR}) we can claim that our proposed polar coding scheme achieves the ergodic capacity of binary-input block fading channels with only CDI.

\subsection{Only CSI-R Available}
From the discussion in Section \ref{S:problemstate} we know that, to achieve capacity in the CSI-R case, it is sufficient to use a polar code for each sub-channel independently with the method of \cite{liu2016fading}. The benefit of independent coding is that different sub-channels can be decoded in parallel, which can improve the throughput of the communication system. Nevertheless, a joint encoding and decoding approach can significantly improve the finite length performance.

Consider a BICM-based approach that the transmitter uses an interleaver to randomize symbols in frame. The equivalent channel seen by the receiver after deinterleaving is an i.i.d. fading channel with known CSI-R if the interleaver is perfect. Thus, the scheme of \cite{liu2016fading} is applicable. If the coherent time is $T_c$ and the number of blocks in a frame is $N$, the total code length will be $T_cN$. Thus, to use standard polar codes, $T_cN$ should be the power of 2

The frame error probability of the parallel scheme can be upper bounded by $$P_e^{Para}\leq T_cO(2^{-N^\beta}).$$
As a comparison, the frame error probability of the BICM-base scheme (assuming a perfect interleaver) satisfies $$P_e^{BICM}\leq O(2^{-(T_cN)^\beta}).$$ Obviously the finite length performance is significantly improved by interleaving, especially when $T_c$ is large. However, the decoding latency of the parallel scheme, assuming $T_c$ decoders working simultaneously, is only $1/(T_c(1+\log T_c/\log N))$ of that for the BICM-based scheme with the same frame size, since the decoding complexity of polar codes is $O(N\log N)$. This shows a trade-off between performance and throughput when used in practice. Since the scheme of \cite{liu2016fading} is capacity-achieving, we can claim that our proposed schemes in this subsection are also capacity-achieving.

\subsection{Full CSI Available}
When both the transmitter and the receiver have perfect CSI, we can simply use standard polar codes for each coherent interval if the coherent time is sufficiently large, as \cite{Santos2013rayleigh} has shown. Otherwise, we can use the CSI-R scheme introduced in the last subsection if the transmitter does not do power control, or design a CSI-R scheme for the equivalent channel after power allocation if it does.

\section{Numerical Results}
\label{S:Numerical}
In the CDI case, the average mutual information (per symbol) of a block fading channel is
\begin{equation}
\label{MI_CDI}
\frac{1}{T_c}I(X_{1:T_c};Y_{1:T_c})=\frac{1}{T_c}\sum_{\mathbf{x}}p(\mathbf{x})\int_{\mathbf{y}} p(\mathbf{y}|\mathbf{x})\log\frac{p(\mathbf{y}|\mathbf{x})}{p(\mathbf{y})}\mathrm{d}\mathbf{y},
\end{equation}
and the mutual information of the $j$th ($j\in[T_c]$) sub-channel is
\begin{equation}
\begin{aligned}
&~~~~I(X_j;Y_{1:T_c}|X_{1:j-1})\\
&=\sum_{x_{1:j}}p(x_{1:{j}})\int_{\mathbf{y}} p(\mathbf{y}|x_{1:{j}})\log\frac{p(\mathbf{y}|x_{1:{j}})}{p(\mathbf{y}|x_{1:{j-1}})}\mathrm{d}\mathbf{y}.
\end{aligned}
\end{equation}

As an example, we assume $h$ follows the Rayleigh distribution with PDF
\begin{equation}
f(h)=\frac{h}{\sigma_h^2}e^{-\frac{h^2}{2\sigma_h^2}}.
\end{equation}
Then the transition PDF of (\ref{transp}) is given by
\begin{equation}
p(\mathbf{y}|\mathbf{x})=\int_{0}^{\infty} \big{(}\prod_{j=1}^{T_c}\frac{1}{\sqrt{2\pi\sigma^2}}e^{-\frac{(y_j-hx_j)^2}{2\sigma^2}}\big{)}\frac{h}{\sigma_h^2}e^{-\frac{h^2}{2\sigma_h^2}}\mathrm{d}h.
\end{equation}

\begin{figure}[tb]
	\centering
	\includegraphics[width=8.5cm]{Tc5.jpg}
	\caption{Achievable rates of binary-input block Rayleigh fading channels.} \label{fig:tc4}
\end{figure}
Since this channel is symmetric, we will only consider uniformly distributed channel inputs. Fig. \ref{fig:tc4} shows a comparison of achievable rates of binary-input AWGN channel, binary-input Rayleigh fading channel with CSI-R, and binary-input block Rayleigh fading channels of different coherent time with only CDI. We can see that as the coherent time increases, the achievable rate with only CDI gets closer and closer to that with CSI-R. It has been shown for several cases (e.g., \cite{peleg1998incoherent,liang2004noncoherent}) that the noncoherent capacity of a block fading channel will approach the coherent capacity as $T_c\rightarrow\infty$. In the considered SNR region (-0.6 db to 0 db), a binary-input block Rayleigh fading channel with 5 symbols per coherent interval has a performance gain about 0.25 db over the i.i.d. fading channel under the CDI assumption, and a performance loss about 0.5 db compared with the CSI-R curve.

\begin{figure}[tb]
	\centering
	\includegraphics[width=8.5cm]{sub5.jpg}
	\caption{Capacity of sub-channels of a block Rayleigh fading channel with coherent time $T_c=5$.} \label{fig:sub4}
\end{figure}
Fig. \ref{fig:sub4} shows the achievable rates of five sub-channels of a block Rayleigh fading channel with coherent time $T_c=5$, compared with that of the original channel and the CSI-R rate. We can see that the achievable rate of a sub-channel increases with its index. This can be intuitively  explained as follows. After decoding a sub-channel, the decoder gains more knowledge about the CSI (although not explicitly shown), and the achievable rates of the following sub-channels become larger. Thus, our proposed CDI scheme can be seen as a realization of the joint channel estimation and data transmission paradigm in noncoherent communications.

\section{Discussion}
\label{S:discuss}
In this paper, we take a coded modulation approach to solve the problem of coding for block fading channels. By viewing transmitted symbols in a coherent block as a supersymbol, we design a multilevel polar coding scheme when only CDI is available, and a parallel scheme as well as a BICM-based scheme when instantaneous CSI is available at the receiver or at both sides. All of these schemes are capacity-achieving.

It is known that in coded modulation, BICM can not achieve coded modulation capacity, because it treats all bit levels equally, which in fact may have some dependency. Similarly, in block fading channels, if we design a BICM-based scheme for the CDI case, it will fail to achieve the capacity since it erases the relation between channel uses in the same coherent block. However, if the receiver knows the instantaneous CSI, each channel use will be indeed independent from others from the receiver's perspective, and the BICM-based scheme can also achieve capacity.

Another approach to deal with the problem of coding for block fading channels is the multiple access channel (MAC) approach, which views a block fading channel of coherent time $T_c$ as a $T_c$-user MAC. We will briefly explain the connection between this approach and ours. From the MAC perspective, code design will be based on MAC polar codes (e.g., \cite{arikan2012sw}). Our proposed MLC-based scheme can be seen as a special case of the MAC-based scheme, i.e., it is equivalent to a MAC polar code designed to achieve a corner point of the achievable rate region of the $T_c$-user MAC. By using other permutations for MAC polarization, one can allocate rates for different "users" more flexibly. However, although this connection is valid, the MAC approach looks unnecessarily complicated for our problem. In the end, it's chain rule, not MAC, that matters for our problem.

Although we have only considered binary inputs in this paper, the proposed schemes can be readily applied to non-binary cases. The idea of this work may also be extended to fading channels with memory. We will leave these to our future work.

\bibliographystyle{IEEEtran}


\section{Introduction}
The study of the mathematical properties of multiloop Feynman integrals
has received increasing attention in the last years both by
the physics and the mathematics communities.
The high precision reached by the experimental measurements
carried out at the LHC, in fact, requires on the theory side the calculation of
multi- (typically two- or three-) loop corrections to various complicated $2 \to 1$,
$2 \to 2$ and $2 \to 3$ processes, including the exchange of
massless and massive virtual particles.
Indeed, unravelling their mathematical structures will be crucial
to handle the complexity of such calculations and, on the long run, it could help
gain a deeper understanding of quantum field theory itself.

Impressive progress has been achieved in the last years in both directions.
The development of new techniques for the calculation
of multiloop Feynman integrals has rendered many previously out-of-reach
calculations feasible. Among these, a fundamental role has been played
by integration-by-parts reductions~\cite{Tkachov:1981wb,Chetyrkin:1981qh,Laporta:2001dd}
and the differential equation method~\cite{Kotikov:1990kg,Remiddi:1997ny,Gehrmann:1999as},
augmented more recently by the use of a canonical basis~\cite{Kotikov:2010gf,Henn:2013pwa}.
\textsl{A posteriori}, a reason for the success of these techniques can be traced back to
the almost concurrent ``discovery'' that the $\epsilon$ expansion of a rather large class of multiloop
(mainly \textsl{massless}) Feynman integrals can be expressed in terms of a very well understood
class of special functions, dubbed
multiple polylogarithms~\cite{Goncharov,Goncharov:1998kja,Remiddi:1999ew,Gehrmann:2000zt}.
In particular the understanding of their algebraic and analytical properties,
also thanks to the development
of the symbol and coproduct formalism~\cite{Goncharov:2005sla,Goncharov:2010jf,Brown:2011ik,Duhr:2011zq,Duhr:2012fh}, together
with stable routines for their evaluation on the whole complex
plane~\cite{Gehrmann:2001pz,Gehrmann:2001jv,Vollinga:2004sn},
have played a crucial role in the last years.
Nevertheless, as it has been known for a long time,
starting at the two-loop order the $\epsilon$ expansion of Feynman integrals can involve new mathematical structures
which lie beyond the realm of multiple polylogarithms. The first and best studied example is that
of the two-loop massive Sunrise graph~\cite{Sabry,Broadhurst:1987ei,Bauberger:1994by,Bauberger:1994hx,Caffo:1998du,Laporta:2004rb,Bloch:2013tra,Remiddi:2013joa,
Adams:2013nia,Adams:2014vja,Adams:2015gva,Adams:2015ydq,Bloch:2016izu,Remiddi:2016gno}.
Its imaginary part in $d=4$ space-time dimensions is known to be expressible as linear combination of
complete elliptic integrals of first and second kind. The real part of the graph can be reconstructed by a dispersive relation and involves integrals over elliptic integrals.
In the last years quite some effort has been spent on studying the new functions
appearing in the Sunrise graph, trying to extend the properties
of multiple polylogarithms to embrace also elliptic generalizations of the latter~\cite{Bloch:2013tra,Adams:2013nia,Adams:2014vja,Adams:2015gva,Adams:2015ydq,Bloch:2016izu}.
Much progress has been made but a complete understanding of these functions is still missing,
in particular regarding their analytic continuation and numerical evaluation.
Moreover, it is not clear whether they can be easily used to express also graphs with
more complicated kinematics (i.e.\ three- and four-point functions).
Last but not least, the interplay of these functions with the simpler multiple polylogarithms
and their iterative structure remains unclear in context with the method of differential equations.
There, simpler integrals appear  as inhomogeneous
terms in the differential equations of more complicated integrals. The solution of these equations
requires on the one hand the solution of the homogeneous equations and, on the other, the integration over
the inhomogeneous piece. Multiple polylogarithms
(and in general Chen iterated integrals~\cite{Chen:1977oja})
are particularly well suited for this scope, as they can be defined as iterative integrals over a set of
differential forms.
 Of course, it is a priori unclear whether such a feature can be expected to
hold for arbitrarily complicated Feynman graphs.

More recently, a new picture has started to emerged, where a generalization of this approach
to more complicated cases becomes possible.
When dealing with a system of coupled differential equations, the most non-trivial
step consists in solving the homogeneous part of the system. In~\cite{Primo:2016ebd} it was shown
that, for arbitrarily complicated cases, integral representations of the homogeneous
solutions can be found by computing the maximal cut of the graph under consideration.
Once the homogeneous solution is known, an integral representation for the
inhomogeneous solution is provided by Euler's variation of constants. Even if usually such integrals
cannot be expressed in terms of known special functions,
from a practical point of view we are particularly interested in obtaining results that allow
for fast and reliable numerical evaluations over the physical phase space.
In this context, it was shown in~\cite{Remiddi:2016gno}
that the study of the imaginary part and of the corresponding
 dispersion relations of Feynman integrals within the differential equations framework
 can facilitate obtaining compact one-fold
integral representations for the two-loop massive Sunrise and the Kite integral\footnote{Note that more recently it was shown that the Kite integral can also be re-written in terms of elliptic polylogarithms~\cite{Adams:2016xah}.}.
There are indications that more complicated Feynman integrals with different and unrelated kinematics
can be casted in a similar form, see for example~\cite{Aglietti:2007as,CaronHuot:2012ab} and more recently a set of two-loop planar double boxes relevant for $H+$jet production~\cite{Bonciani:2016qxi}.

The new ideas summarized above provide us with the tools required to start successfully
studying more examples
of relevant Feynman diagrams whose $\epsilon$ expansions do not evaluate to multiple polylogarithms.
In this paper, we focus on a three-point non-planar topology
which is relevant for the two-loop QCD
 corrections to the processes $gg\to t\bar{t}$
 and $g g \to \gamma \gamma$
 through a massive top loop.
Similar integrals
appear in the non-planar sectors of $H+$jet and $HH$ production.
The rest of the paper is organized as follows.
In Section~\ref{sec:notation} we describe the problem and establish the notation.
We continue in Section~\ref{sec:subtopos} showing how the $\epsilon$ expansion
of the simpler subtopologies can be integrated in terms of multiple polylogarithm
with a new algorithm based on the differential equations in canonical form.
In Section~\ref{sec:toptopo} we consider the top topology, which cannot be integrated in terms
of multiple polylogarithms. We show how to solve the coupled system of
differential equations satisfied by its two master integrals based on
information extracted from the maximal cut;
we then compare it to the result obtained by solving directly the second order differential equation.
We thoroughly study the required functions in Section~\ref{sec:toptopo}
and Appendix~\ref{sec:ancont}, where
we also show how to analytically continue our solution to all relevant regions of the phase space.
Finally in Section~\ref{sec:inhomsol} we use the results of the previous sections
to build up the inhomogeneous solution
for the top topology in analytic form. We write the results
as one-fold integrals over rational and irrational functions,
multiple polylogarithms and complete elliptic integrals of first and second kind.
We then draw our conclusions in Section~\ref{sec:conclusions}.

\section{Differential equations}\label{sec:notation}
We consider the Feynman integrals family
$I_{a_1,a_2,a_3,a_4,a_5,a_6,a_7}$ defined as
\begin{align}
& \trianglecross{q}{p_1}{p_2} = I_{a_1,a_2,a_3,a_4,a_5,a_6,a_7} \Big|_{a_7 < 0} \nonumber \\
& = \int \,
\frac{ \mathrm{d}^dk\, \mathrm{d}^d l \;\; (m^2)^{2\epsilon} N_\epsilon^2 \;\; (k^2)^{-a_7}}{[(k-p_1)^2]^{a_1} [(l-p_1)^2-m^2]^{a_2} [(k+p_2)^2]^{a_3}
[(k-l+p_2)-m^2]^{a_4} [(k-l)-m^2]^{a_5}[(l^2-m^2)^{a_6}]}\,, \label{eq:topo}
\end{align}
where
\begin{equation}
N_\epsilon = \frac{1}{i \pi^{2-\epsilon} \Gamma(1+\epsilon)}
\end{equation}
and $\epsilon=(4-d)/2$.
We stress that the triangle topology above is obtained by limiting ourselves to negative
powers of the last propagator, i.e.\ $a_7<0$.
Four of the six propagators are massive (thick lines), while only one of the three external legs is
off-shell, i.e.\ $p_1^2 = p_2^2 = 0$ and $q^2 =(p_1+p_2)^2 = s$\,.
Integrals in family~\eqref{eq:topo} with different values of the exponents $a_i$
may be reduced to master integrals using integration by parts reductions
as implemented for example in Reduze\;2~\cite{vonManteuffel:2012np,Studerus:2009ye,Bauer:2000cp,fermat}.
As a result we find 11 different master integrals: 9 for the subtopologies
and 2 for the top topology.

As we will see explicitly, the nine master integrals of the subtopologies
can be expressed in terms of multiple polylogarithms only.
For them we choose a canonical basis
using the method described in~\cite{Gehrmann:2014bfa}.\footnote{For recent related work on the analysis of master
integrals and the construction of a canonical basis see
\cite{Lee:2013hzt,Argeri:2014qva,Lee:2014ioa,Tancredi:2015pta,Ita:2015tya,
Meyer:2016slj,Georgoudis:2016wff,Prausa:2017ltv,Gituliar:2017vzm}.
}
The remaining two integrals, instead, contain elliptic generalizations of the multiple polylogarithms and, therefore, it is not possible to find for them a canonical basis as defined in~\cite{Henn:2013pwa}.\footnote{We stress here that
a canonical basis as defined in~\cite{Henn:2013pwa} requires not only a specific factorization of
the $\epsilon$ dependence from the kinematics, but also, as an additional condition on the matrix of the system,
a d-$\log$ form for the differential equations.
In the general case the first condition could still be realized, while the d-$\log$ form
condition might have to be generalized to differentials of more complicated functions.}
This can be motivated using the ideas presented in~\cite{Tancredi:2015pta}. There it was shown
for various examples that the possibility of decoupling the differential equations in the limit $d \to 4$ (and therefore of
writing the result in terms of multiple polylogarithms) is signaled by a degeneracy
of the integration by parts identities in the limit of even numbers of dimensions,
$d \to 2\,n$ with $n \in \mathbb{N}$.

We first consider the part of the Euclidean region where $s<-4\,m^2$ and employ the
following basis
\begin{align}
m_{1} &= \epsilon^2\, I^{2,10}_{0,2,0,2,0,0,0}\,, &
m_{2} &= \epsilon^2\, s\, I^{3,7}_{2,2,1,0,0,0,0} \nonumber \\
m_{3} &= \epsilon^2\,s\, I^{3,22}_{0,2,1,0,2,0,0}\,, &
m_{4} &= \epsilon^2\, \sqrt{s(s-4\,m^2)}\, \left[ I^{3,22}_{0,2,2,0,1,0,0}
            + \frac{1}{2}  I^{3,22}_{0,2,1,0,2,0,0} \right] \nonumber \\
m_{5} &= \epsilon^3\, s\, I^{4,15}_{1,2,1,1,0,0,0}\,, &
m_{6} &= \epsilon^2\,\sqrt{s(s+4\,m^2)} \left[ s \, I^{4,15}_{2,2,1,1,0,0,0}
            - \frac{\epsilon}{2\,m^2 (1+2 \epsilon)}\, I^{2,10}_{0,2,0,2,0,0,0}\right] \nonumber \\
m_{7} &= \epsilon^3\, s\, I^{4,30}_{0,2,1,1,1,0,0}\,, &
m_{8} &= \epsilon^4\, s\, I^{5,59}_{1,1,0,1,1,1,0}\,, \qquad\;
m_{9} = \epsilon^4 \, s\, I^{5,31}_{1,1,1,1,1,0,0}\,, \nonumber \\
m_{10} &= \epsilon^4 s^2 I^{6,63}_{1, 1, 1, 1, 1, 1, 0}\,, &
m_{11} &= \epsilon^4 \frac{s^2 (s + 16 m^2)}{2 (1 + 2 \epsilon)} I^{6,63}_{1, 2, 1, 1, 1, 1, 0}\,,
\label{eq:basis}
\end{align}
where we added the number of different denominators and {Reduze}'s sector id as a superscript to
the $I$ integrals for easier identification.

In order to write down the system of differential equations we introduce the
massless ratio
\[
x = -s/m^2\,.
\]
We introduce the vector $\vec{m} = (m_1, \ldots, m_9)$ for the master
integrals of the subtopologies, which fulfil canonical differential
equations. In matrix notation we obtain
\begin{equation} \label{eq:deqcan}
\mathrm{d} \vec{m} = \epsilon\, \sum_{i=1}^5 \mathrm{d}\,\ln(l_i(x))\,A_i \,\vec{m}
\end{equation}
with the letters
\begin{equation}\label{eq:alphabet}
 l_1 = \sqrt{x},\quad
 l_2 = \tfrac{1}{2}(\sqrt{x} + \sqrt{x+4}),\quad
 l_3 = \sqrt{x + 4},\quad
 l_4 = \tfrac{1}{2}(\sqrt{x} + \sqrt{x - 4}),\quad
 l_5 = \sqrt{x - 4}
\end{equation}
and the matrices
\begin{alignat}{3}
A_1 &= \left( \begin{smallmatrix}
0& 0& 0& 0& 0& 0& 0              & 0& 0\\
0&-2& 0& 0& 0& 0& 0              & 0& 0\\
0& 0& 2& 0& 0& 0& 0              & 0& 0\\
0& 0& 0&-2& 0& 0& 0              & 0& 0\\
0& 0& 0& 0& 2& 0& 0              & 0& 0\\
0& 0& 0& 0& 0&-2& 0              & 0& 0\\
0& 0& 2& 0& 0& 0& 0              & 0& 0\\
0& 0& 0& 0& 0& 0& 2              & 0& 0\\
0& 0& 0& 0& 2& 0&-1& 0& -2
\end{smallmatrix} \right), &\quad
A_2 &= \left( \begin{smallmatrix}
0& 0& 0& 0& 0& 0& 0& 0& 0\\
0& 0& 0& 0& 0& 0& 0& 0& 0\\
0& 0& 0& 4& 0& 0& 0& 0& 0\\
2& 0& -6& 0& 0& 0& 0& 0& 0\\
0& 0& 0& 0& 0& 0& 0& 0& 0\\
0& 0& 0& 0& 0& 0& 0& 0& 0\\
0& 0& 0& 0& 0& 0& 0& 0& 0\\
0& 0& 0& 0& 0& 0& 0& 0& 0\\
0& 0& 0& 0& 0& 0& 0& 0& 0
\end{smallmatrix} \right), &\quad
A_3 &= \left( \begin{smallmatrix}
0& 0& 0& 0& 0& 0& 0& 0& 0\\
0& 0& 0& 0& 0& 0& 0& 0& 0\\
0& 0& 0& 0& 0& 0& 0& 0& 0\\
0& 0& 0&-6& 0& 0& 0& 0& 0\\
0& 0& 0& 0& 0& 0& 0& 0& 0\\
0& 0& 0& 0& 0& 0& 0& 0& 0\\
0& 0& 0& 0& 0& 0& 0& 0& 0\\
0& 0& 0& 0& 0& 0& 0& 0& 0\\
0& 0& 0& 0& 0& 0& 0& 0& 0
\end{smallmatrix} \right),
\nonumber\\
A_4 &= \left( \begin{smallmatrix}
0& 0& 0& 0& 0& 0& 0& 0& 0\\
0& 0& 0& 0& 0& 0& 0& 0& 0\\
0& 0& 0& 0& 0& 0& 0& 0& 0\\
0& 0& 0& 0& 0& 0& 0& 0& 0\\
0& 0& 0& 0& 0& 2& 0& 0& 0\\
2& 2& 0& 0&-6& 0& 0& 0& 0\\
0& 0& 0& 0& 0& 0& 0& 0& 0\\
0& 0& 0& 0& 0& 0& 0& 0& 0\\
0& 0& 0& 0& 0& 0& 0& 0& 0
\end{smallmatrix} \right), &\quad
A_5 &= \left( \begin{smallmatrix}
0& 0& 0& 0& 0& 0& 0& 0& 0\\
0& 0& 0& 0& 0& 0& 0& 0& 0\\
0& 0& 0& 0& 0& 0& 0& 0& 0\\
0& 0& 0& 0& 0& 0& 0& 0& 0\\
0& 0& 0& 0& 0& 0& 0& 0& 0\\
0& 0& 0& 0& 0&-4& 0& 0& 0\\
0& 0& 0& 0& 0& 0& 0& 0& 0\\
0& 0& 0& 0& 0& 0& 0& 0& 0\\
0& 0& 0& 0& 0& 0& 0& 0& 0
\end{smallmatrix} \right).
&&
\end{alignat}
The root appearing in the normalization of $m_4$ leads to $l_2$ and $l_3$,
while the root appearing in the normalization of $m_6$ leads to $l_4$ and $l_5$.
We note that, as long as we limit ourselves to $m_1,\ldots,m_9$,
these pairs of roots never mix and we could rationalize them separately
with two different changes of variables.
The two structures, nevertheless, mix up once considering the differential equations for the
two master integrals of the top sector, $m_{10}$ and $m_{11}$.

With the basis given in~\eqref{eq:basis}, the differential equations for the master integrals
of the top topology read
\begin{align}
\frac{d}{dx} \left(\begin{matrix} m_{10} \\ m_{11} \end{matrix} \right) =
B(x) \left(\begin{matrix} m_{10} \\ m_{11} \end{matrix} \right)
+ \epsilon\, D(x) \left(\begin{matrix} m_{10} \\ m_{11} \end{matrix} \right)
+ \left(\begin{matrix} N_{10}(\epsilon;x) \\ N_{11}(\epsilon;x) \end{matrix} \right) \label{eq:topsys}
\end{align}
where $B(x)$ and $D(x)$ are two $2 \times 2$ matrices that do not depend on $\epsilon$,
while $N_{10}(\epsilon;x)$ and $N_{11}(\epsilon;x)$  contain the dependence on the
simpler subtopologies.
The matrices of the homogeneous part are
\begin{align}
B(x) = \left( \begin{array}{ccc} 0 & & \frac{1}{2(x-16)} - \frac{1}{2\,x} \\&\\
                                        \frac{1}{2\,x} & & \frac{1}{x}  \end{array}\right)\,,
                                        \qquad
D(x) = \left( \begin{array}{ccc} - \frac{2}{x} & & \frac{1}{x-16} - \frac{1}{x} \\&\\
                                        \frac{2}{x} & & \frac{1}{x} - \frac{1}{x-16} \end{array}\right)\,,
\end{align}
while the non-homogeneous terms read
\begin{align}
&N_{10}(\epsilon; x) = 0\,, \nonumber \\
&N_{11}(\epsilon; x) = \frac{\epsilon^2}{1+2\epsilon} \left( \frac{5}{4} \,m_3 - \frac{7}{2} \,m_5
- \frac{3}{2} \frac{x}{\sqrt{x(x-4)}} \, m_6 + 2\, m_7 + m_9\, \right)\,. \label{eq:nonhom}
\end{align}
From the latter it is clear that in the top topology all letters
appear simultaneously.
In our normalization the integrals $m_j = m_j(\epsilon; x)$ $(j=1,\ldots,11)$ have a
Taylor series expansion at $\epsilon=0$,
\begin{equation}\label{mepsexpansion}
m_j(\epsilon; x) = \sum_{n=0}^{\infty} m_j^{(n)}(x)\; \epsilon^n\qquad\text{for}~j=1,\ldots,11,
\end{equation}
and we wish to calculate the expansion coefficients up to $n=4$.

\section{Integration of the subtopologies}\label{sec:subtopos}

We first consider the evaluation of the subtopologies $m_1,\ldots,m_9$.
Inserting the expansion~\eqref{mepsexpansion} into the differential
equation~\eqref{eq:deqcan} and comparing coefficients of powers in $\epsilon$ gives
\begin{equation}\label{deqexpansion}
\mathrm{d} m_j^{(n)}(x) = \sum_{i=1}^5 \mathrm{d}\,\ln(l_i(x))\,A_i \,m_j^{(n-1)}(x)\qquad\text{for~}n=0,\ldots,\infty,
\end{equation}
which allows us to solve fully decoupled systems of differential equations
order by order $n$ bottom-up.
The differential equations~\eqref{deqexpansion} have the form of the
differential equations of multiple polylogarithms with root-valued
letters.
We find it useful to solve the full vector $\vec{m}$ in a uniform setup,
where we consider all letters at the same time.
This prevents a rationalization of the $x$ dependence and
therefore a traditional integration of the differential equations.

Instead, we construct an ansatz for the solution using suitable target
functions and require it to fulfil the differential equation~\eqref{deqexpansion}.
In other words, we \emph{integrate the symbol} of the loop integrals.
The basis of our function construction is the Duhr-Gangl-Rhodes algorithm~\cite{Duhr:2011zq},
which we extend here to the case of an irrational alphabet.
For simpler cases it is typically not strictly necessary to work with
a canonical basis in order to be able to integrate the differential equations,
provided they (partially) decouple.
In contrast to this, our strategy here crucially relies on the
differential equation~\eqref{deqexpansion} being in canonical form,
such that we can read off the symbol letters of the solution. A similar approach
for the integration of the differential equations has been used in~\cite{Bonciani:2016qxi}
to calculate solutions through to weight two.

To construct our solutions in the region $s<-4m^2$ through to weight four
we choose the functions
\begin{equation}
\ln,\quad \Li_2,\quad \Li_3,\quad \Li_4,\quad \Li_{2,2}
\end{equation}
with suitable arguments.
We look for functions which do not introduce symbol letters beyond the alphabet
\eqref{eq:alphabet} of the differential equation~\eqref{deqexpansion}.
For the arguments of $\ln$ we can just choose the original letters themselves.

A classical polylogarithm $\Li_j(z)$ has symbol entries $z$ and $1-z$.
Choosing $z$ to be a power product of our original letters,
\begin{equation}
\label{powerproduct}
c^{a_0} l_1^{a_1} \cdots l_5^{a_5},
\end{equation}
with $c$ a rational number, ensures $z$ to factorize over our alphabet.
Here, it is enough to consider $a_n\in \mathbbm{Z}$ and $c^{a_0}=\pm 1$.
We stress that this is too constraining in general.
In particular, depending on the (choice of) letters it may be necessary
to also consider roots of letters for factorization and allow the $a_n$ to be
non-integer rational numbers.
Note that there is an ambiguity related to the freedom of redefining the
letters by expressions
of the form \eqref{powerproduct}.
In our case, the definition of our letters absorbs any further explicit
roots.

In addition to the argument $z$ of a polylogarithm we also require
$1-z$ to factorize in the sense \eqref{powerproduct} over our original alphabet.
This condition is challenging to check with computer algebra systems, in particular
for multivariate applications.
We therefore resort to decompositions found by inserting numbers
for the variables and applying heuristic integer relation searches to the logarithms
of the involved factors~\cite{weihsprivcomm,Lenstra1982,PARI2} (more details will be given in~\cite{rmsavm}).
Our heuristic search produces this admissible set:
\begin{align}
\{&1/(l_1 l_4),\; l_4/l_1,\; 1/(l_2 l_3),\; l_3/
 l_2,\; -(1/(l_1 l_2)),\; -(1/(l_4 l_5)),\; l_5/l_4,\; 1/l_4^2,\; l_4/
 l_5,\; l_2/l_3,\; l_2/l_1,\; l_4^2,\; 1/
 l_2^2,\;  \nonumber\\
 &-(1/l_4^2),\;
-l_4 l_5,\; -(1/l_2^2),\; l_2 l_3,\; l_1/l_4,\;
 l_1 l_4,\; -l_4^2,\; l_1/l_2,\;
 l_2^2,\; -l_2^2,\; -l_1 l_2,\; -(1/(l_1 l_2^2 l_3)),\; -(1/(
  l_1 l_4^2 l_5)),\;\nonumber\\
 & l_4^2/(l_1 l_5),\; 1/l_2^4,\;
 l_2^2/(
 l_1 l_3),\; 1/l_4^4,\; l_4^4,\; (l_1 l_5)/l_4^2,\; -l_1 l_4^2 l_5,\;
 l_2^4,\; (l_1 l_3)/l_2^2,\; -l_1 l_2^2 l_3\},
\end{align}
which is closed under $z \to 1-z$ and $z \to 1/z$, where $z$ is an element of the set.
For $\Li_{2,2}(z_1,z_2)$ we need to make sure that the symbol factors
$z_1$, $1-z_1$, $z_2$, $1-z_2$ and $z_1-z_2$ do not introduce new letters.
Our construction proceeds along the lines described above, but we
do not list the generic result for the argument set here since it is too lengthy.

We concentrate here on the part of the Euclidean region where $s<-4m^2$
and require all functions to be real valued for $x>4$.
Indeed, we find that all of these constraints can be satisfied and a solution
in terms of such functions can be found, which satisfies the differential
equations.
Our functions read
\begin{align}
\{
 &\ln(l_1), \ln(l_2), \ln(l_3), \ln(l_4), \ln(l_5),
 \Li_2(l_2^{-2}), \Li_2(1/(l_2 l_3))), \Li_2(-l_4^{-2}), \Li_2(-1/(l_4 l_5)),
 \Li_3(l_2^{-2}),
 \nonumber \\
 &\Li_3(l_1/l_2), \Li_3(-1/(l_1 l_2^2 l_3)),
 \Li_3(1/(l_2 l_3))), \Li_3(l_2/l_3)), \Li_3(-l_4^{-2}), \Li_3(l_4^{-2}),
 \Li_3(1/(l_1 l_4)),
 \nonumber \\
 & \Li_3(-1/(l_4 l_5))), \Li_3(l_1 l_5)/l_4^2),
 \Li_4(-l_2^{-2}), \Li_4(l_2^{-2}), \Li_4(-1/(l_1 l_2)), \Li_4(l_1/l_2),
 \Li_4(-1/(l_1 l_2^2 l_3)),
 \nonumber \\
 &
 \Li_4(1/(l_2 l_3)), \Li_4(l_2/l_3), \Li_4(l_1 l_3)/l_2^2), \Li_4(-l_4^{-2}), \Li_4(l_4^{-2}),
 \Li_4(1/(l_1 l_4)), \Li_4(l_4/l_1), \Li_4(-1/(l_1 l_4^2 l_5)),
 \nonumber \\
 &
 \Li_4(-1/(l_4 l_5)), \Li_4(l_1 l_5/l_4^2), \Li_4(l_5/l_4), \Li_{2,2}(-1, -l_2^{-2}),
 \Li_{2,2}(-1, -l_4^{-2}), \Li_{2,2}(-1/(l_1 l_2), l_1/l_2),
 \nonumber \\
 &
 \Li_{2,2}(1/(l_1 l_4), l_1/l_4)
\}\,.
\end{align}
In order to fix the integration constants, we use the massive tadpole
and the massless bubbles as input as well as various regularity conditions.

We find very compact results for the final solutions of the subtopology
integrals in terms of our choice of functions.
These results will be a necessary building block to construct the full solution of the top topology.
The explicit form of the results can be found in Appendix~\ref{app:subtopos}
and, in electronic format, in the ancillary file of the arXiv submission of this paper.

\section{Integration of the top topology}\label{sec:toptopo}

\subsection{Homogenous solution}

Having the expressions for all subtopologies in terms of multiple polylogarithms,
we can now proceed and study the differential equations~\eqref{eq:topsys} for the top sector.
We start by expanding in $\epsilon$ the inhomogeneous terms~\eqref{eq:nonhom}, up to
$\epsilon^4$
\begin{equation}
N_{j}(\epsilon;x) = \sum_{a=0}^4 \, N_j^{(a)}(x)\,\epsilon^a\,,\qquad \mbox{for}\;j=10,11\,.
\end{equation}
Substituting explicitly the subtopologies we find the remarkably simple expressions
\begin{align}
&N_{10}^{(0)}(x) = 0\,, \quad N_{10}^{(1)}(x) = 0\,, \quad N_{10}^{(2)}(x) = 0\,,\quad N_{10}^{(3)}(x) =0\, \quad
N_{10}^{(4)}(x) = 0\,,\nonumber \\
&N_{11}^{(0)}(x) = 0\,, \quad N_{11}^{(1)}(x) = 0\,, \quad N_{11}^{(2)}(x) = 0\,,\quad N_{11}^{(3)}(x) =0 \,,\nonumber \\
&N_{11}^{(4)}(x) = 5 \ln^2 (l_2)
- l_1\; \frac{3/2\, \zeta_2
+ 3 \ln^2( l_4) + 3 \Li_2( -1/l_4^2 )}{ l_5}\,, \label{eq:nonhomexp}
\end{align}
where the letters $l_j$ have been defined above~\eqref{eq:alphabet}.
We can now expand the right- and left-hand-side of~\eqref{eq:topsys} and we find,
at any order $n$,
\begin{align}
\frac{d}{dx} \left(\begin{matrix} m^{(n)}_{10} \\ m_{11}^{(n)} \end{matrix} \right) =
B(x) \left(\begin{matrix} m^{(n)}_{10} \\ m^{(n)}_{11} \end{matrix} \right)
+  D(x) \left(\begin{matrix} m^{(n-1)}_{10} \\ m^{(n-1)}_{11} \end{matrix} \right)
+ \left(\begin{matrix} N^{(n)}_{10}(x) \\ N^{(n)}_{11}(x) \end{matrix} \right)\,.\label{eq:topsysex}
\end{align}
Comparing with~\eqref{eq:nonhomexp}, we see that we get a non trivial
non-homogeneous term only at order $\epsilon^4$, which is also the maximum
order in the expansion we are interested in.
In fact, the integrals $I^{6,63}_{1, 1, 1, 1, 1, 1, 0}$
and $I^{6,63}_{1, 2, 1, 1, 1, 1, 0}$ are finite, which can be checked
e.g.\ with the methods of \cite{Panzer:2014gra,vonManteuffel:2014qoa} implemented in {Reduze\;2.1},
such that the first non-vanishing coefficients of $m_{10}$ and $m_{11}$ occur at order $\epsilon^4$.

The first step to solve~\eqref{eq:topsysex} is to find a solution of its homogeneous part,
which at every order in $\epsilon$ is a coupled system of two first-order linear differential equations
\begin{equation}
\frac{d}{dx} \left(\begin{matrix} M_{10} \\ M_{11} \end{matrix} \right) =
B(x) \left(\begin{matrix} M_{10} \\ M_{11} \end{matrix} \right)\,, \label{eq:tophom}
\end{equation}
where we use a capital letter $M$ instead of $m$ to describe the homogeneous part
of the solution and drop the superscript $(n)$, since the form of the equation does not depend on $n$.
In order to solve the system we need to find two independent set of solutions, i.e.\ a $2 \times 2$ matrix
\begin{align}
G(x) &= \left(\begin{matrix} I_1(x) & J_1(x)  \\ I_2(x) & J_2(x) \end{matrix}\right) \label{eq:matsol}
\end{align}
such that
\begin{equation}
\frac{d}{dx} G(x) = B(x) \, G(x)\,. \label{eq:deqmat}
\end{equation}
The inverse of the matrix $G(x)$ reads
\begin{equation}
G^{(-1)}(x) = \frac{1}{W(x)} \left( \begin{array}{cc} J_2(x) & -J_1(x) \\ -I_2(x) & I_1(x)  \end{array}\right)\,, \label{eq:invG}
\end{equation}
where $W(x)$ is the Wronskian of the solutions, $W(x) = I_1(x)J_2(x) - J_1(x)I_2(x)$\,.
If the matrix $G(x)$ is known, one can perform at every order in $\epsilon$ the rotation
\begin{align}
\left(\begin{matrix} m^{(n)}_{10} \\ m^{(n)}_{11} \end{matrix} \right)
&= G \left(\begin{matrix} \tilde{m}^{(n)}_{10} \\ \tilde{m}^{(n)}_{11} \end{matrix} \right) \label{eq:rotation}
\end{align}
such that the new functions $\tilde{m}^{(n)}_{10}, \tilde{m}^{(n)}_{11}$ fulfil the differential equations
\begin{align}
\frac{d}{dx} \left(\begin{matrix} \tilde{m}^{(n)}_{10} \\ \tilde{m}_{11}^{(n)} \end{matrix} \right) =
G^{(-1)}(x)\,D(x)\, G(x) \left(\begin{matrix} \tilde{m}^{(n-1)}_{10} \\ \tilde{m}^{(n-1)}_{11} \end{matrix} \right)
+ G^{(-1)}(x) \left(\begin{matrix} N^{(n)}_{10}(x) \\ N^{(n)}_{11}(x) \end{matrix} \right)\, \label{eq:toprot}
\end{align}
whose solution can be now written by quadrature as
\begin{align}
\left(\begin{matrix} \tilde{m}^{(n)}_{10}(y) \\ \tilde{m}_{11}^{(n)}(y) \end{matrix} \right) =
\int_{y_0}^y \, dx\, G^{(-1)}(x)\,D(x)\, G(x) \left(\begin{matrix} \tilde{m}^{(n-1)}_{10} \\ \tilde{m}^{(n-1)}_{11} \end{matrix} \right)
+ \int_{y_0}^y \, dx\, G^{(-1)}(x) \left(\begin{matrix} N^{(n)}_{10}(x) \\ N^{(n)}_{11}(x) \end{matrix} \right)
+ \left(\begin{matrix} c^{(n)}_{10} \\ c^{(n)}_{11} \end{matrix} \right) \label{eq:soltoprot}
\end{align}
where $c_{10}^{(n)}$ and $c_{11}^{(n)}$ are two integration constants and $y_0$ is a suitable boundary for the integration.

In order to find an explicit solution for the matrix $G(x)$ it is useful to recast~\eqref{eq:tophom} as a
second order differential equation for the first master integral
\begin{equation}
\frac{d^2 \, M_{10}(x) }{dx^2} + \left( \frac{1}{x-16}\right)\frac{d \, M_{10}(x) }{dx}
+ \frac{1}{64}\left( \frac{1}{x}  + \frac{16}{x^2} - \frac{1}{x-16} \right)M_{10}(x) = 0\,. \label{eq:secdeqm}
\end{equation}
A general way to solve this equation is to compute the maximal cut of $M_{10}(x)$
as explained in detail in~\cite{Primo:2016ebd}. In that reference it was shown that
the maximal cut of $M_{10}(x)$ computed on a specific integration contour can be written
(up to an irrelevant overall numerical constant) as
\begin{align}
{\rm Cut_1}(M_{10}(x)) = \sqrt{ \frac{x}{x-16} } \EK\left( \frac{x}{x-16}\right)\,, \label{eq:cut1}
\end{align}
where we introduced the complete elliptic integral of the first kind defined as
\begin{equation}
\EK(z) = \int_0^1 \frac{dx}{\sqrt{(1-x^2)(1-\,z\,x^2)}} \,, \quad
\mbox{for} \quad z \in \mathbb{C} \quad \mbox{and}\quad \Re{(z)} < 1\,. \label{eq:ellint}
\end{equation}
It is straightforward to verify that~\eqref{eq:cut1} satisfies~\eqref{eq:secdeqm}.
A second independent solution can be obtained integrating on an independent contour.
As explained in~\cite{Primo:2016ebd}, this can be avoided in the case of elliptic integrals,
as the second solution is simply given by
\begin{align}
{\rm Cut_2}(M_{10}(x)) = \sqrt{ \frac{x}{x-16} } \EK\left( 1- \frac{x}{x-16}\right)\,. \label{eq:cut2}
\end{align}
The two solutions~\eqref{eq:cut1} and \eqref{eq:cut2} completely determine the first row of the matrix~\eqref{eq:matsol}.

Before proceeding, it is instructive to try to solve Eq.~\eqref{eq:secdeqm} directly, since in this case the equation is particularly
simple. We define the auxiliary function
\begin{equation}
f(x) = \frac{1}{\sqrt{x}}\, M_{10}(x)\,,
\end{equation}
which fulfills the differential equation
\begin{equation}
\frac{d^2 \, f(x) }{dx^2} + \left( \frac{1}{x} + \frac{1}{x-16} \right) \frac{d \,f (x) }{dx}
+ \frac{1}{64} \left(\frac{1}{x-16} - \frac{1}{x} \right)\, f(x) = 0 \,. \label{eq:secdeq}
\end{equation}
Equation~\eqref{eq:secdeq} is now in standard form and its solution can be
easily written in terms
of the complete elliptic integral of the first kind~\eqref{eq:ellint}.
Since the differential equation~\eqref{eq:secdeq} has regular singular points in $x=0$, $x=16$
and $x=\pm \infty$, we must consider the solution of the latter in the three intervals $0<x<16$,
$16<x<\infty$ and $-\infty<x<0$\,.
We start by considering the region $0<x<16 $\, i.e.\ $-16\,m^2 < s < 0$.
This region is non physical and the master integrals $m_{10}$, $m_{11}$ are real there.
Following the same procedure, we can build up similar solutions in the remaining
regions, i.e.\ the part of the Euclidean region where  $-\infty < s < -16\,m^2$ and
the physical region $s>0$.
We derive explicitly those solutions in Appendix~\ref{sec:ancont}.
We use the notation $I_j^{(a,b)}(x)$ and $J_j^{(a,b)}(x)$, with $j=1,2$, for the two pairs
of solutions valid in the region $a<x<b$, and consequently $G^{(a,b)}(x)$ for the matrix of the solutions
in the same region.
The general solution of the homogeneous second order differential equation for $0<x<16$
can be written as
\begin{equation}
\label{eq:fsol}
f(x) = c_{1}\, \EK \left(  \frac{x}{16} \right) + c_{2} \EK \left( 1 - \frac{x}{16} \right)\,,
\end{equation}
which determines the first row of the matrix~\eqref{eq:matsol} as
\begin{align}
I_1^{(0,16)}(x) = \sqrt{x}\, \EK \left(  \frac{x}{16} \right)\,, \qquad J_1^{(0,16)}(x) =
\sqrt{x}\, \EK \left( 1 - \frac{x}{16} \right)\,. \label{eq:sol1r1}
\end{align}

Eqs.~\eqref{eq:sol1r1} are apparently different from the solutions determined by studying
the maximal cut~\eqref{eq:cut1} and \eqref{eq:cut2}. Indeed, as the four functions are all solutions
of the same second order homogeneous differential equation, a relation among them must
exist such that only two functions are really independent.
Indeed, it is well known that the elliptic integral of the first kind $\EK(z)$ and its complementary
$\EK(1-z)$ fulfil the following relations
\begin{align}
&\EK\left( z\right) = \frac{1}{\sqrt{z}} \left( \EK\left( \frac{1}{z}\right) -   i\, \EK\left( 1-\frac{1}{z}\right) \right)  \label{eq:relell1} \\
&\EK\left( 1-z\right) = \frac{1}{\sqrt{z}}  \EK\left( 1-\frac{1}{z}\right)
\label{eq:relell2}
\end{align}
where we use the prescription
$z \to z + i\, 0^+$ for definiteness.\footnote{
Note that Eqs.~\eqref{eq:relell1} and \eqref{eq:relell2} are not
the only relations between elliptic integrals of different arguments.
Other relations are known, which we did
not find necessary in this context.}
Applying Eq.~\eqref{eq:relell1} twice to the solutions~\eqref{eq:cut1} and \eqref{eq:cut2},
we construct two more equivalent sets of solutions
\begin{align}
\left\{ \EK\left(\frac{16}{x}\right)\,, \;\;  \EK\left( 1- \frac{16}{x}\right) \right\} \qquad
\mbox{and} \qquad
\left\{ \sqrt{x}\,  \EK\left(\frac{x}{16}\right)\, ,\; \; \sqrt{x}\,\EK\left( 1- \frac{x}{16}\right) \right\}.
\label{eq:invariantset}
\end{align}
It is simple to check by direct insertion that these functions do indeed solve
our second order differential equation.
This proves the equivalence of the solutions found through the maximal cut with the ones found
by solving directly the differential equation.
The existence of more than one representation for the same solutions in terms of elliptic integrals
of different arguments will turn out to be very important to write down the analytic continuation
of the solution from $0<x<16$ to the whole phase space in terms of explicitly real solutions,
as explained in Appendix~\ref{sec:ancont}.

We can now proceed with the solution of the system~\eqref{eq:topsysex}
in the region $0<x<16$; here the most convenient representation is provided
by the choice~\eqref{eq:sol1r1}, which we will adopt from now on.
In order to determine the second row of the matrix $G(x)$ we
plug~\eqref{eq:sol1r1} into~\eqref{eq:deqmat} and find for consistency that
\begin{align}
I_2^{(0,16)}(x) = -\sqrt{x}\,\EE \left(  \frac{x}{16}  \right)\,, \qquad J_2^{(0,16)}(x) =
\sqrt{x} \left[ \EE \left( 1-\frac{x}{16} \right) - \EK \left( 1- \frac{x}{16} \right) \right]\,, \label{eq:sol2r1}
\end{align}
where we introduced the complete elliptic integral of the second kind
\begin{equation}
\EE(z) = \int_0^1 dx\, \frac{\sqrt{1-\,z\,x^2}}{\sqrt{1-x^2}} \,, \quad
\mbox{for} \quad z \in \mathbb{C} \quad \mbox{and}\quad \Re{(z)} < 1\,. \label{eq:ellint2}
\end{equation}
We write therefore the matrix of the solutions valid for $0<x<16$ as
\begin{align}
G^{(0,16)}(x) = \left(\begin{matrix} I_1^{(0,16)}(x) & J_1^{(0,16)}(x)  \\ I_2^{(0,16)}(x) & J_2^{(0,16)}(x) \end{matrix}\right)\,. \label{eq:solmatr1}
\end{align}
We can now compute the Wronskian of the solutions. We find
\begin{align}
\label{eq:wronskianval}
W^{(0,16)}(x) = \det{ \left( G^{(0,16)}(x) \right) } =
 I_1^{(0,16)}(x) J_2^{(0,16)}(x) - J_1^{(0,16)}(x) I_2^{(0,16)}(x) = \frac{\pi}{2}\,x
 \end{align}
which is a direct consequence of the Legendre identity among the first two complete elliptic integrals
\begin{equation}
\EK(z) \EE(1-z) + \EE(z) \EK(1-z) - \EK(z) \EK(1-z) = \frac{\pi}{2}\,.
\end{equation}

Alternatively, we can determine $W(x)$ also without resorting to the Legendre identity.
The fact that the Wronskian must be a linear function of $x$ can be easily seen taking its derivative as
\begin{equation}
\frac{d}{dx} W(x) = \frac{d}{dx} \det{\left( G(x) \right)} =
Tr\left( G^{(-1)}(x) B(x) G(x)\right) \det \left( G(x) \right) = \frac{1}{x}  W(x)\,,
\end{equation}
which gives as a solution
\begin{equation}
W(x) = c\, x\,, \label{eq:linwr}
\end{equation}
with $c$ an arbitrary integration constant.\footnote{Note that this is of course a general property
of the solution and remains true
independent of the region $a<x<b$ that we consider.}
The value of the constant $c$ can be then fixed
by computing the Wronskian for a fixed value of $x$.
For example, we can study the behaviour of the solutions at the two boundaries, i.e.\
$x \to 0^+$ and $x \to 16^-$. We find respectively, keeping only the leading behaviour,

\begin{align}
\lim_{x \to 0^+} \; I_1^{(0,16)}(x) &= \frac{\pi}{2} \sqrt{x} \,, \nonumber\\
\lim_{x \to 0^+} \; I_2^{(0,16)}(x) &= - \frac{\pi}{2} \sqrt{x} \,,  \nonumber \\
\lim_{x \to 0^+} \; J_1^{(0,16)}(x) &= \frac{\sqrt{x}}{2}\left( -\ln{(x)} + 8\, \ln{(2)} \right)\,, \nonumber \\
\lim_{x \to 0^+} \; J_2^{(0,16)}(x) &= \frac{\sqrt{x}}{2}\left( \ln{(x)} - 8\, \ln{(2)} + 2 \right) \,,\label{eq:expr1a}
\end{align}
and
\begin{align}
\lim_{x \to 16^-} \; I_1^{(0,16)}(x) &= -2\ln{(16-x)} + 16 \ln{(2)}
+ \frac{(16-x)}{32} \left( \ln{(16-x)} - 8 \ln{(2)} - 2\right) \,,\nonumber \\
\lim_{x \to 16^-} \; I_2^{(0,16)}(x) &=  -4
+ \frac{(16-x)}{16} \left( \ln{(16-x)} - 8 \ln{(2)} + 3\right)\,,\nonumber \\
\lim_{x \to 16^-} \; J_1^{(0,16)}(x) &= 2\,\pi  - \frac{\pi}{32} (16-x) \,,\nonumber \\
\lim_{x \to 16^-} \; J_2^{(0,16)}(x) &= - \frac{\pi}{16} (16-x) \,. \label{eq:expr1b}
\end{align}
By using~\eqref{eq:expr1a} and~\eqref{eq:expr1b} it is easy to verify explicitly that
\begin{align}
\lim_{x \to 0^+} \frac{W^{(0,16)}(x)}{x} = \lim_{x \to 16^-} \frac{W^{(0,16)}(x)}{x} = \frac{\pi}{2}\,,
\end{align}
which fixes the constant $c$.

\subsection{Analytic continuation of the inhomogeneous term}

With the results of the previous sections we can solve the inhomogeneous
differential equations for the top topology using Euler's
variation of constants.
In order to employ explicitly real building blocks, we need to analytically
continue the homogeneous solutions and the inhomogeneous term to the various
regions of the phase space.
Details for the continuation of the homogeneous solutions are worked out
in Appendix~\ref{sec:ancont}.
For the analytic continuation of the inhomogeneous term we give the explicit
results in this section to highlight the different emerging functions.

For clarity, in every region $a<x<b$ we
extract explicitly the imaginary part (whenever there is one) and
introduce the functions $\mathcal{R}^{(a,b)}(x)$ and  $\mathcal{Q}^{(a,b)}(x)$
respectively for the real and imaginary parts of $N_{11}^{(4)}(x)$ in that region
$$N_{11}^{(4)}(x) \Big|_{a<x<b} = \mathcal{R}^{(a,b)}(x) + i\, \mathcal{Q}^{(a,b)}(x)\,.$$
For the analytical continuation to other kinematical regions we assign an
infinitesimal positive imaginary part to $s$ which translates to a negative
imaginary part for $x$, $x \to x - i\,0^+$.
We will also make use of letters from the alphabet
\begin{align}
l'_i &=  \{ \sqrt{-x}, \tfrac{1}{2}(\sqrt{-x} + \sqrt{-x-4}), \sqrt{-x - 4},
\tfrac{1}{2}(\sqrt{-x} + \sqrt{-x + 4}), \sqrt{-x + 4} \},
\end{align}
depending on the region of phase space.

For the region $4<x<16$ we obtain
\begin{align}
 \mathcal{R}^{(4,16)}(x) &= 5 \ln^2 (l_2)
- 3\, l_1\; \frac{ \zeta_2/2
+  \ln^2( l_4) +  \Li_2( -1/l_4^2 )}{ l_5}\,,
\nonumber\\
\mathcal{Q}^{(4,16)}(x) &= 0\,. \label{eq:R416}
\end{align}
We remark that $x=16$ is actually no special point for the subtopologies
and a single analytic expression is sufficient for the entire range $x>4$.
In the region $0<x<4$, $N_{11}^{(4)}(x)$ stays real and we obtain
\begin{align}
\mathcal{R}^{(0,4)}(x) &= 5 \ln^2 (l_2)
+ 3\,l_1\; \frac{\Cl_2(-2 \arccsc(2/l_1) )}{ l_5'}\,,
\nonumber\\
\mathcal{Q}^{(0,4)}(x) &=0\,, \label{eq:R04}
\end{align}
where the Clausen function
$
\Cl_2(\theta) = - \int_0^\theta \ln{\left| 2 \sin{\frac{t}{2}}\right|}\, dt
$
arises from the dilogarithm of a pure phase factor,
$\Cl_2(\theta) = \Im\left(\Li_2(e^{i \theta})\right)$,
which can be seen from $l_4^2 = (\sqrt{x}-i\sqrt{4-x})/(\sqrt{x}+i\sqrt{4-x})$.
In the region $-4 < x < 0$, the inhomogeneous terms develops an imaginary part and can be
conveniently expressed in terms of
\begin{align}
\mathcal{R}^{(-4,0)}(x) &= - 5\arcsec^2{(2/l_3)} + 3\, l_1'\, \frac{\left(   \zeta_2  -\ln^2{(l_4')}   -
\Li_2(1/l_4'^2) \right)}{l_5'}\,,
\nonumber\\
\mathcal{Q}^{(-4,0)}(x) &= \pi\, \frac{3\, l_1'}{l_5'} \ln{(l_4')}\,, \label{eq:Rm40}
\end{align}
where the first term of the real part arises from a purely imaginary $\ln(l_2)$
since $l_2^2 = (\sqrt{-x} + i \sqrt{4 + x})/(\sqrt{-x} - i \sqrt{4 + x})$
is a pure phase factor.
Finally, for $-\infty < x < -4$ we obtain for the real and the imaginary part,
\begin{align}
\mathcal{R}^{(-\infty,-4)}(x) &= 5\ln^2{(l_2')} - \frac{15}{2}\zeta_2 + 3\, l_1'\, \frac{\left(  \zeta_2 - \ln^2{(l_4')}   - \, \Li_2(1/l_4'^2) \right)}{l_5'}\,,
\nonumber\\
\mathcal{Q}^{(-\infty,-4)}(x) &= \pi \,\left( \frac{3\, l_1'}{l_5'} \ln{(l_4')}- 5 \, \ln{(l_2')} \right)\,, \label{eq:Rmoom4}
\end{align}
respectively.

\subsection{The inhomogeneous solution}\label{sec:inhomsol}

We are now finally in the position of writing down the complete solution
of~\eqref{eq:topsys} in all relevant regions of the phase space.
For definiteness, we start in the Euclidean region $0<x<4$, $-4 m^2 < s < 0$,
and then use the results of the previous section and of Appendix~\ref{sec:ancont} to continue the solution
to the other regions.
In Appendix~\ref{sec:ancont}, in particular, we build up matrices of solutions
 $$G^{(a,b)}(x) =
 \left( \begin{matrix} I_1^{(a,b)}(x) && J_1^{(a,b)}(x) \\
 I_2^{(a,b)}(x) && J_2^{(a,b)}(x) \end{matrix}\right)$$
whose entries are explicitly real for $a<x<b$, and we show explicitly how to match
the homogeneous solution moving from one region to the other.
We recall that the expansions of $m_{10}$ and $m_{11}$ start at order $\epsilon^4$.
We first write down the solutions for $\tilde{m}_{10}^{(4)}$ and $\tilde{m}_{11}^{(4)}$ and then rotate
them back to obtain expressions for the physical master integrals $m_{10}$ and $m_{11}$.

\subsubsection*{Region $\mathbf{0<x<16}$}

From~\eqref{eq:toprot}, \eqref{eq:invG}, \eqref{eq:wronskianval}
and the inhomogeneous term in~\eqref{eq:nonhomexp}
we get in the region $0<x<4$
\begin{align}
\frac{d}{dx}\,\left( \begin{matrix} \tilde{m}_{10}^{(4)} \\ \tilde{m}_{11}^{(4)} \end{matrix}\right) =
\frac{2}{\pi} \frac{1}{x}\,\left( \begin{matrix} J_2^{(0,16)}(x) && -J_1^{(0,16)}(x) \\ -I_2^{(0,16)}(x) && I_1^{(0,16)}(x) \end{matrix}\right)
\left( \begin{matrix} 0 \\  \mathcal{R}^{(0,4)}(x) \end{matrix}\right)
\end{align}
which gives
\begin{align}
&\tilde{m}_{10}^{(4)}(x) = - \frac{2}{\pi}\,\int_{0}^{x} \, \frac{dy}{y}\, J_1^{(0,16)}(y)\, \mathcal{R}^{(0,4)}(y) + c_{10}
\,, \nonumber\\
&\tilde{m}_{11}^{(4)}(x) =  + \frac{2}{\pi}\, \int_{0}^x \, \frac{dy}{y}\, I_1^{(0,16)}(y)\,  \mathcal{R}^{(0,4)}(y) + c_{11}\,,
\label{eq:solmtilde}
\end{align}
where $c_{10}$ and $c_{11}$ are two integration constants.
We stress again that $\mathcal{R}^{(0,4)}(y)$ is a simple combination of weight two
polylogarithms.
Since neither the top topology nor the subtopologies develop
an imaginary part when $x$ crosses the pseudo-threshold $x=4$,
we can write an explicitly real solution for the physical (unrotated) master integrals
$m_{10}$, $m_{11}$
valid for the larger range $0<x<16$
\begin{align}
m_{10}^{(4)}(x) &=  c_{10}\, I_1^{(0,16)}(x) + c_{11}\, J_1^{(0,16)}(x) \nonumber \\
&+\frac{2}{\pi} \, J_1^{(0,16)}(x)\, \int_{0}^x \,\frac{dy}{y}\, I_1^{(0,16)}(y)\,\mathcal{R}^{(a,b)}(y)
- \frac{2}{\pi} \, I_1^{(0,16)}(x)\, \int_{0}^x \, \frac{dy}{y}\, J_1^{(0,16)}(y)\, \mathcal{R}^{(a,b)}(y)
\,, \nonumber\\ & \nonumber \\
m_{11}^{(4)}(x) &=  c_{10}\, I_2^{(0,16)}(x) + c_{11}\, J_2^{(0,16)}(x)\,, \nonumber \\
&+\frac{2}{\pi} \, J_2^{(0,16)}(x)\, \int_{0}^x \, \frac{dy}{y}\, I_1^{(0,16)}(y)\, \mathcal{R}^{(a,b)}(y)
- \frac{2}{\pi} \, I_2^{(0,16)}(x)\, \int_{0}^x \, \frac{dy}{y}\, J_1^{(0,16)}(y)\, \mathcal{R}^{(a,b)}(y)\,,
\label{eq:soltofix}
\end{align}
where $(a,b)$ is either $(0,4)$ or $(4,16)$, depending on whether $y<4$ or $y>4$.
The integration constants must be fixed with two proper boundary conditions.
Imposing that the integrals be regular as $x \to 16^-$
and go to zero for $x \to 0^+$, we find $c_{10}=c_{11}=0$
\begin{align}
m_{10}^{(4)}(x) &=  \frac{2}{\pi} \, \left[ J_1^{(0,16)}(x)\, \int_{0}^x \, \frac{dy}{y}\, I_1^{(0,16)}(y)\,
\mathcal{R}^{(a,b)}(y)
- I_1^{(0,16)}(x)\, \int_{0}^x \, \frac{dy}{y}\, J_1^{(0,16)}(y)\, \mathcal{R}^{(a,b)}(y)
\right]
\,, \nonumber\\ & \nonumber \\
m_{11}^{(4)}(x) &=
\frac{2}{\pi} \, \left[ J_2^{(0,16)}(x)\, \int_{0}^x \, \frac{dy}{y}\, I_1^{(0,16)}(y)\, \mathcal{R}^{(a,b)}(y)
- I_2^{(0,16)}(x)\, \int_{0}^x \, \frac{dy}{y}\, J_1^{(0,16)}(y)\, \mathcal{R}^{(a,b)}(y) \right] \,.
\label{eq:solfixed}
\end{align}

\subsubsection*{Region $\mathbf{-4<x<0}$}

We continue the solution to the region $-4<x<0$ using $x \to x - i\,0^+$,
\eqref{eq:Rm40}, \eqref{eq:matching} and~\eqref{eq:mxmatch} and obtain
\begin{align}
m_{10}^{(4)}(x) =   \frac{2}{\pi} \,
&\left[ J_1^{(-\infty,0)}(x)\, \int_{0}^x \, \frac{dy}{y}\, I_1^{(-\infty,0)}(y)\,  \mathcal{R}^{(-4,0)}(y)
  -  \, I_1^{(-\infty,0)}(x)\, \int_{0}^x \, \frac{dy}{y}\, J_1^{(-\infty,0)}(y)\,  \mathcal{R}^{(-4,0)}(y) \right]\nonumber \\
+ i\, \frac{2}{\pi}\,
&\left[  J_1^{(-\infty,0)}(x)\, \int_{0}^x \, \frac{dy}{y}\, I_1^{(-\infty,0)}(y)\, \mathcal{Q}^{(-4,0)}(y)
-  I_1^{(-\infty,0)}(x)\, \int_{0}^x \, \frac{dy}{y}\, J_1^{(-\infty,0)}(y) \, \, \mathcal{Q}^{(-4,0)}(y)\right]
\,, \nonumber\\  & \nonumber \\
m_{11}^{(4)}(x) =   \frac{2}{\pi} \,
&\left[ J_2^{(-\infty,0)}(x)\, \int_{0}^x \, \frac{dy}{y}\, I_1^{(-\infty,0)}(y)\, \mathcal{R}^{(-4,0)}(y)
-  \, I_2^{(-\infty,0)}(x)\, \int_{0}^x \, \frac{dy}{y}\, J_1^{(-\infty,0)}(y)\, \mathcal{R}^{(-4,0)}(y) \right]\nonumber \\
+ i\, \frac{2}{\pi}\,
&\left[  J_2^{(-\infty,0)}(x)\, \int_{0}^x \, \frac{dy}{y}\, I_1^{(-\infty,0)}(y)\, \mathcal{Q}^{(-4,0)}(y)
-  I_2^{(-\infty,0)}(x)\, \int_{0}^x \, \frac{dy}{y}\, J_1^{(-\infty,0)}(y) \, \, \mathcal{Q}^{(-4,0)}(y)\right]\,.
\label{eq:solmink1}
\end{align}
Note that now all integrals are explicitly real and all imaginary parts are explicit.

\subsubsection*{Region $\mathbf{-\infty<x<-4}$}
We can then go further and continue beyond the threshold for the production of two massive
particles, i.e.\ $s>4 m^2$, $x<-4$. Using again formulas~\eqref{eq:Rmoom4},~\eqref{eq:matching} and~\eqref{eq:mxmatch}
we obtain
\begin{align}
m_{10}^{(4)}(x) =    \frac{2}{\pi} \,
&\left[ J_1^{(-\infty,0)}(x)\, \left( \int_{-4}^x \, \frac{dy}{y}\, I_1^{(-\infty,0)}(y)\,  \mathcal{R}^{(-\infty,-4)}(y) -
 c_1 \right) \right. \nonumber \\
  -  \, &\;\left. I_1^{(-\infty,0)}(x)\, \left( \int_{-4}^x \, \frac{dy}{y}\, J_1^{(-\infty,0)}(y)\,  \mathcal{R}^{(-\infty,-4)}(y)
  -  c_2 \right) \right]\nonumber \\
+ i\, \frac{2}{\pi}\,
&\left[  J_1^{(-\infty,0)}(x)\,\left( \int_{-4}^x \, \frac{dy}{y}\, I_1^{(-\infty,0)}(y)\,  \mathcal{Q}^{(-\infty,-4)}(y) -
 c_3 \right) \right. \nonumber \\
  -  \, &\;\left. I_1^{(-\infty,0)}(x)\, \left( \int_{-4}^x \, \frac{dy}{y}\, J_1^{(-\infty,0)}(y)\,  \mathcal{Q}^{(-\infty,-4)}(y)
  -  c_4\right)\right]
\,,
\nonumber\\
m_{11}^{(4)}(x) =    \frac{2}{\pi} \,
&\left[ J_2^{(-\infty,0)}(x)\, \left( \int_{-4}^x \, \frac{dy}{y}\, I_1^{(-\infty,0)}(y)\,  \mathcal{R}^{(-\infty,-4)}(y) -
 c_1 \right) \right. \nonumber \\
  -  \, &\;\left. I_2^{(-\infty,0)}(x)\, \left( \int_{-4}^x \, \frac{dy}{y}\, J_1^{(-\infty,0)}(y)\,  \mathcal{R}^{(-\infty,-4)}(y)
  -  c_2 \right) \right]\nonumber \\
+ i\, \frac{2}{\pi}\,
&\left[  J_2^{(-\infty,0)}(x)\,\left( \int_{-4}^x \, \frac{dy}{y}\, I_1^{(-\infty,0)}(y)\,  \mathcal{Q}^{(-\infty,-4)}(y) -
 c_3 \right) \right. \nonumber \\
  -  \, &\;\left. I_2^{(-\infty,0)}(x)\, \left( \int_{-4}^x \, \frac{dy}{y}\, J_1^{(-\infty,0)}(y)\,  \mathcal{Q}^{(-\infty,-4)}(y)
  -  c_4\right)\right]
\,,
\label{eq:solmink2}
\end{align}
where all imaginary parts are explicit and all integrals are real.
The matching constants are
\begin{align}
&c_1 = \int_{-4}^0 \, \frac{dy}{y}\, I_1^{(-\infty,0)}(y)\,  \mathcal{R}^{(-4,0)}(y)\,, \quad
c_2 = \int_{-4}^0 \, \frac{dy}{y}\, J_1^{(-\infty,0)}(y)\,  \mathcal{R}^{(-4,0)}(y) \,, \nonumber \\
&c_3 = \int_{-4}^0 \, \frac{dy}{y}\, I_1^{(-\infty,0)}(y)\,  \mathcal{Q}^{(-4,0)}(y) \,, \quad
c_4 =  \int_{-4}^0 \, \frac{dy}{y}\, J_1^{(-\infty,0)}(y)\,  \mathcal{Q}^{(-4,0)}(y) \,.
\end{align}
They can be evaluated with high precision starting
directly from this integral representation.

\subsubsection*{Numerical evaluation}

We verified our formulas in the previous section with a thorough comparison
against SecDec\;3, both in the
physical and in the non-physical region.
Numerical results for the top-level master integrals are shown in Figure~\ref{fig:i1x11110}.
We would like to point out that our representation allows for particularly fast
and precise evaluations also in the physical region of phase space.
As an example for a result with higher precision, we obtain
\begin{align}
\left. m^4 I_{1,1,1,1,1,1,0}\right|_{s=5m^2, d=4} &\approx
-0.07776462028160023644086669458011467822536257409024
\nonumber\\
&~+ i~ 0.34306740464518688969054397597465622650767181505054
\end{align}
for one of our six-line master integrals.
\begin{figure}[ht]
\centering
  \includegraphics[width=\linewidth]{i1x11110}
  \caption{\label{fig:i1x11110}
  Finite top-level master integrals $m^4\, I_{1,1,1,1,1,1,0} $ (left)
  and $m^6\, I_{1,2,1,1,1,1,0} $ (right).
  Solid curves were obtained with our representations,
  dots with SecDec\;3.
  }
\end{figure}

\section{Conclusions} \label{sec:conclusions}

Feynman integrals which evaluate to classes of functions outside the realm of
multiple polylogarithms
constitute the bottleneck for many multiloop calculations relevant for LHC phenomenology.
To this day, only a very limited number of examples have been considered in the literature.
An important step towards a more complete understanding of the new mathematical
structures therefore consists of the study of more explicit examples in order to expose their common properties.
In this paper we have considered the calculation of a two-loop non-planar three-point function
with four massive internal propagators and one external off-shell leg
using the method of differential equations.
The differential equations for the subtopologies can be put in canonical form
and integrated in terms of multiple polylogarithms over an irrational alphabet.
In general, in the presence of different roots, traditional approaches to
integrate the result in terms of Goncharov's polylogarithms fail,
unless one is able to find a change of variables which linearizes
or at least rationalizes the letters of the alphabet.
We employed a new algorithm to perform the integration without any
rationalization of the alphabet.
This allowed us to write the results up to weight four in terms
of simple logarithms, classical polylogarithms and $\Li_{2,2}$ functions only.

We moved then to the two master integrals of the top level sector, which fulfil
two coupled differential equations in the momentum transfer squared. By studying the maximal cut
of the two master integrals we showed how to solve the differential equations
in terms of complete elliptic integrals of first and second kind; we then used Euler's variation
of constants in order to write the solution as a one-fold integral over elliptic integrals and
multiple polylogarithms.
Finally, we used the same techniques described in~\cite{Remiddi:2016gno} in order to
continue our results to the Minkowski physical region. In this way we provided
expressions which can be evaluated in a fast and reliable manner in the
entire physical region of the phase space.
We compared our results with {SecDec\;3} and found perfect agreement.
We expect that similar techniques can be extended to study also more complicated
three- and four-point functions, as in part confirmed by the results
recently obtained in~\cite{Bonciani:2016qxi}.

\section*{Acknowledgments}

We would like to thank Stefan Groote, Ettore Remiddi, Robert Schabinger
and Erich Weihs for valuable discussions.
We are grateful to the Mainz Institute for Theoretical Physics (MITP) for its
hospitality and support.