\section{Introduction}
\label{sec:Intro}

In mathematical chemistry as well as in graph theory there is a great interest in the study of various graph parameters that are defined as sums of certain quantities over all vertices (such as the first Zagreb index and the Graovac-Pisanski index), over all pairs of adjacent vertices (such as the Szeged index, the second Zagreb index, the Balaban index, the Randi\'c index, the Mostar index and the irregularity of a graph) or over all pairs of vertices of the corresponding graph (such as the Wiener index and the total irregularity of a graph). See, for example, \cite{BDHKS, AD, DMSTZ18, GFK, KST, KST1, WB} for the most recent research on these indices. In many of these parameters the above mentioned quantities are distance-related. For instance, in the Wiener index the corresponding quantity is the distance between the two vertices in question while in the Szeged index the corresponding quantity for the pair $\{u,v\}$ is the product $n_{u,v}n_{v,u}$, where for any two vertices $x$ and $y$ of the graph in question $n_{x,y}$ is the number of vertices of the graph that are closer to $x$ than to $y$.

The most recently introduced index of the above mentioned ones is the {\em Mostar index}~\cite{DMSTZ18} which is defined as
$$
	\mathrm{Mo}(\Gamma) = \sum_{\{u,v\} \in E(\Gamma)}|n_{u,v} - n_{v,u}|.
$$
This index is intimately related to the concept of distance-balancedness of graphs, which was first studied in~\cite{Han99} although the term {\em distance-balanced} was introduced nine years later~\cite{JKR08}. In terms of the above mentioned Mostar index a graph is distance-balanced if and only if its Mostar index is equal to $0$.

In~\cite{MikSpa18} the following related concept was introduced. For a pair of (not necessarily adjacent) vertices $u$ and $v$ of a (connected) graph we say that this pair of vertices is {\em balanced} if $n_{u,v} = n_{v,u}$. Thus, a connected graph is distance-balanced if and only if every pair of adjacent vertices is balanced. As there is no good reason why one would only be interested in the question of when every pair of adjacent vertices in a graph is balanced, the authors of~\cite{MikSpa18} (but see also~\cite{FreMik18}) proposed to study the following generalizations of the notion of distance-balancedness. A graph of diameter at least $\ell$ is said to be {\em $\ell$-distance-balanced} if each pair of its vertices at distance $\ell$ is balanced. Therefore, $1$-distance-balancedness is equivalent to being distance-balanced in the sense of~\cite{JKR08}.

One could thus generalize the Mostar index in the sense that the {\em $\ell$-th Mostar index} of a connected finite graph $\Gamma$ would be defined as
$$
	\mathrm{Mo}^\ell(\Gamma) = \sum_{\{u,v\} \subseteq V(\Gamma) \atop d(u,v) = \ell}|n_{u,v} - n_{v,u}|.
$$
However, to measure how much the graph is (un)balanced (in terms of distances) as a whole, one should really study the following parameter:
$$
	\mathrm{Mo}^*(\Gamma) = \mathrm{uB}(\Gamma) = \sum_{\{u,v\} \subseteq V(\Gamma)} |n_{u,v} - n_{v,u}| = \sum_{\ell = 1}^D \mathrm{Mo}^\ell(\Gamma),
$$
where $D$ is the diameter of $\Gamma$. We call $\mathrm{Mo}^*(\Gamma) = \mathrm{uB}(\Gamma)$ the {\em distance-unbalancedness (index)} of $\Gamma$. In~\cite{MikSpa18} the notion of high distance-balancedness was introduced where a graph is {\em highly distance-balanced} if each pair of vertices of the graph is balanced. In other words, $\Gamma$ is highly distance-balanced if and only if $\mathrm{uB}(\Gamma) = 0$. Therefore, the index $\mathrm{uB}(\Gamma)$ is a measure of how ``unbalanced'' in terms of distances the graph $\Gamma$ is.

As with the Wiener index (which is defined as the sum of distances between all pairs of vertices of the graph), where this number divided by the number of all distinct pairs of vertices of the graph is also studied, it might also be interesting to study the {\em average distance-unbalancedness} of a graph. For a graph $\Gamma$ of order $n$ this number is defined as
$$
	\overline{\mathrm{uB}}(\Gamma) = \frac{\mathrm{uB}(\Gamma)}{{n \choose 2}}.
$$

There are of course numerous interesting natural questions that arise. How close to $0$ can the (average) distance-unbalancedness of a graph that is not highly distance-balanced be? For each positive integer $n$ what is the smallest nonzero and the largest possible value of $\mathrm{uB}(\Gamma)$ when $\Gamma$ runs through the set of all connected graphs of order $n$? What if $\Gamma$ runs just through the set of all members of order $n$ of a certain family of graphs (such as trees)?

Our paper is organized as follows. After some preliminaries in Section \ref{sec:Prelim} we compute the distance-unbalancedness index of some well-known families in Section \ref{sec:Fam}. We consider distance-unbalancedness of trees and general graphs in Sections \ref{sec:trees} and \ref{sec:min}, respectively. Along the way, a few conjectures are stated and some open problems are proposed.

\section{Preliminaries}
\label{sec:Prelim}

Throughout this paper all graphs are simple, undirected, finite and connected (unless otherwise specified). We first introduce some additional notation that will shorten the arguments regarding balancedness of pairs of vertices. For a graph $\Gamma$ and its vertices $u, v$ we let $W_{u,v}(\Gamma)$ be the set of all vertices of $\Gamma$ that are closer to $u$ than to $v$, and we let $n_{u,v}(\Gamma)$ be its cardinality (if the graph $\Gamma$ is clear from the context we simply write $W_{u,v}$ and $n_{u,v}$). The {\em distance-unbalancedness} of the pair $\{u,v\}$ is then defined as $\mathrm{uB}_\Gamma(u,v) = |n_{u,v}(\Gamma) - n_{v,u}(\Gamma)|$ (again, we usually omit the subscript $\Gamma$). The Mostar index of a connected graph $\Gamma$ is thus simply $\mathrm{Mo}(\Gamma) = \sum_{\{u,v\} \in E(\Gamma)}\mathrm{uB}(u,v)$ and the distance-unbalancedness index of $\Gamma$ is $\mathrm{uB}(\Gamma) = \sum_{\{u,v\} \subseteq V(\Gamma)}\mathrm{uB}(u,v)$.

For a positive integer $n$ the complete bipartite graph $K_{n,1}$ will be called the {\em star} of order $n+1$ and will be denoted by $S_n$. The graph obtained by taking a disjoint union of stars $S_n$ and $S_m$ and then adding an edge between the central vertex of $S_n$ (the one of valence $n$) and the central vertex of $S_m$ will be denoted by $\SS{n}{m}$ and called a {\em merged star}. The graph obtained from $\SS{n}{m}$ by subdividing the edge between the central vertices of $S_n$ and $S_m$ with one vertex (usually denoted by $x$) will be denoted by $\SSx{n}{m}$ and called a {\em subdivided merged star}.

The cycle (graph) of order $n \geq 3$ will be denoted by $C_n$ and the path of order $n$ (and thus length $n-1$) will be denoted by $P_n$.

For an integer $k \geq 3$ and positive integers $n_1 \geq n_2 \geq \cdots \geq n_k$ we let the {\em spider graph} $\mathrm{Sp}(n_1,n_2,\ldots , n_k)$ be the tree of order $1+n_1+n_2+\cdots + n_k$ having one vertex of valence $k$ and all other vertices of valence at most $2$ such that the graph obtained by removing this unique vertex of valence $k$ is a disjoint union of $k$ paths of orders $n_1, n_2, \ldots , n_k$. For instance, for $n \geq 3$ the star $S_n$ is in fact the spider graph $\mathrm{Sp}(1,1,\ldots, 1)$ of order $n+1$.

\section{Distance-unbalancedness of some well-known families}
\label{sec:Fam}

As was pointed out in~\cite{MikSpa18} a sufficient condition for a pair of vertices $\{u,v\}$ of a graph $\Gamma$ to be balanced is that there exists an automorphism of $\Gamma$ interchanging $u$ and $v$. This simple but useful fact will be used throughout the paper. As an immediate corollary (see~\cite[Corollary~3.4]{MikSpa18}) each Cayley graph of an Abelian group is highly distance-balanced and thus its distance-unbalancedness index is $0$.

In what follows we compute the distance-unbalancedness index $\mathrm{uB}(\Gamma)$ for members of some well-known families of graphs. We first consider complete multipartite graphs.

\begin{proposition}
\label{pro:multi}
Let $n_1 \geq n_2 \geq n_3 \geq \cdots \geq n_k$ be positive integers. Then
$$
	\mathrm{uB}(K_{n_1, n_2, \ldots , n_k}) = \sum_{1 \leq i < j \leq k}n_in_j(n_i - n_j).
$$
\end{proposition}

\begin{proof}
Let $\Gamma = K_{n_1, n_2, \ldots , n_k}$. Partition its vertex-set into sets $V_1, V_2, \ldots , V_k$ with $|V_i| = n_i$, $1 \leq i  \leq k$, in such a way that for each $i$, $1 \leq i  \leq k$, each vertex of $V_i$ is adjacent to every vertex outside $V_i$ (and to no vertex from $V_i$). As any pair of vertices from $V_i$, $1 \leq i \leq k$, can be interchanged by an automorphism of $\Gamma$, the only possible unbalanced pairs of vertices are pairs $\{u,v\}$, where $u \in V_i$ and $v \in V_j$ for some $1 \leq i < j \leq k$. Clearly, $n_{u,v} = |1+(n_j-1) - (1 + (n_i-1))| = n_i - n_j$, and so the result follows.
\end{proof}

As an immediate corollary of the above proposition we get the distance-unbalancedness index of stars. Namely, for any positive integer $n \geq 1$ we get that $\mathrm{uB}(S_n) = \mathrm{uB}(K_{n,1}) = n(n-1)$ and consequently $\overline{\mathrm{uB}}(S_n) = \frac{2(n-1)}{n+1}$. The {\em wheel} graphs, which are related to them (the wheel graph $W_{n}$, where $n \geq 3$, is obtained from the star $S_n$ by connecting all of its leafs into a single cycle) have a slightly smaller distance-unbalancedness index. Namely, the following holds (the proof of the stated result is straightforward and is left to the reader).

\begin{proposition}
For each $n \geq 3$ we have that $\mathrm{uB}(W_n) = n(n-3)$.
\end{proposition}

Even though we will have a closer look at distance-unbalancedness of trees in the next section we now compute this index for a very special kind of trees.

\begin{proposition}
\label{pro:paths}
Let $n$ be a positive integer. Then the following holds for the path $P_{n}$ of order $n$:
$$
	\mathrm{uB}(P_{n}) = \left\{\begin{array}{cc}
		\frac{(n-2)n(2n+1)}{12} & \text{if }  n\ \text{is even},\\
		\frac{(n-1)(n+1)(2n-3)}{12} &\text{if }  n\ \text{is odd},\end{array}\right. \quad \text{and}\quad
		\overline{\mathrm{uB}}(\Gamma) = \left\{\begin{array}{ccc}
		\frac{(n-2)(2n+1)}{6(n-1)} &\text{if }  n\ \text{is even},\\
		\frac{(n+1)(2n-3)}{6n} &\text{if }  n\ \text{is odd}.\end{array}\right.
$$
\end{proposition}

\begin{proof}
Denote the vertices of $\Gamma = P_{n}$ by $v_1, \ldots, v_n$, where $v_i \sim v_{i+1}$ for all $i$, $1 \leq i < n$. Let $i,j$ be integers with $1 \leq i < j \leq n$ and let $\Gamma'$ be the subgraph of $\Gamma$ induced on the set of all vertices $v_k$ with $i \leq k \leq j$. As there is an automorphism of $\Gamma'$ interchanging $v_i$ and $v_j$, it follows that $\mathrm{uB}_{\Gamma'}(v_i, v_j) = 0$, and so it clearly follows that $\mathrm{uB}_\Gamma(v_i,v_j) = |(i-1)-(n-j)| = |n+1-i-j|$. Therefore,
$$
	\mathrm{uB}(\Gamma) = \sum_{1 \leq i < j \leq n}|n+1-i-j| = \sum_{i = 1}^{n-1}\sum_{j=i+1}^n |n+1-i-j|.
$$
Let us denote this number by $a_n$. To compute $a_n$ one thus has to sum up all the numbers from the following upper triangular table (where the lines correspond to all possible values of $i$ from $1$ to $n-1$ in increasing order):
$$
	\begin{array}{ccccccccc}
	n-2 & n-3 & n-4 & n-5 & \cdots & 3 & 2 & 1 & 0 \\
	     & n-4 & n-5 & n-6 & \cdots & 2 & 1 & 0 & 1 \\
	     &      & n-6 & n-7 & \cdots & 1 & 0 & 1 & 2 \\
	     &      &      &      &  \vdots &   &   &   &  \\
	     &      &      &      &           & n-8 & n-7  & n-6  & n-5 \\
     	     &      &      &      &           &      & n-6  & n-5  & n-4 \\
     	     &      &      &      &           &      &       & n-4  & n-3 \\
     	     &      &      &      &           &      &       &       & n-2 \end{array}
$$
It is thus clear that $a_{n+2} = a_n + n(n+1)$ holds for all $n \geq 0$ where of course $a_1 = a_2 = 0$. Solving this elementary recursion yields the result of this proposition.
\end{proof}

We can use the above result to calculate the distance-unbalancedness index of some tube-like graphs (Cartesian products of paths by cycles). Note that each Cartesian product of two cycles is highly distance-balanced as it is a Cayley graph of an abelian group (the direct product of the two corresponding cyclic groups) and thus has distance-unbalancedness index $0$. This of course also holds for trivial tubes $P_1 \square C_m$ and $P_2 \square C_m$.

We first fix some notation pertaining to tube-like graphs. Let $n \geq 1$ and $m \geq 3$ be integers and let $\Gamma = P_n \square C_m$. The vertices of $\Gamma$ will be denoted by pairs $(i,j)$, $i \in \{1,2,\ldots , n\}$, $j \in \mathbb{Z}_m$ in the natural way (so that for each $1 \leq i \leq n$ and each $j \in \mathbb{Z}_m$ the vertex $(i,j)$ is adjacent to $(i,j\pm 1)$, and for each $1 < i < n$ and each $j \in \mathbb{Z}_m$ the vertex $(i,j)$ is adjacent to $(i\pm 1, j)$). For each $i \in \{1,2,\ldots , n\}$ we let $\Gamma^i$ be the subgraph of $\Gamma$ induced on the set of vertices $\{(i,j) : j \in \mathbb{Z}_m\}$. We call $\Gamma^i$ the {\em $i$-th column} of $\Gamma$. Similarly, for each $j \in \mathbb{Z}_m$ the $j$-th {\em row} $\Gamma_j$ of $\Gamma$ is the subgraph of $\Gamma$ induced on the set of vertices $\{(i,j) \colon 1 \leq i \leq n\}$.

Since the permutations $\alpha, \beta, \gamma$, mapping each vertex $(i,j)$ to $(i,j+1)$, $(i,-j)$ and $(n-i+1,j)$, respectively (where the second coordinate is computed modulo $m$), are automorphisms of $\Gamma$, it is clear that for any $j, j' \in \mathbb{Z}_m$ there is an automorphism of $\Gamma$ interchanging the vertices $(1,j)$ and $(n,j')$, and so the corresponding pair of vertices is balanced in $\Gamma$. A similar argument shows the following. Let $1 \leq i_1 < i_2 \leq n$ and let $j_1, j_2 \in \mathbb{Z}_m$. Then $(i_1,j_1)$ and $(i_2,j_2)$ are balanced in the subgraph of $\Gamma$, induced by the union of columns $\Gamma^{i_1} \cup \Gamma^{i_1+1} \cup \cdots \cup \Gamma^{i_2}$. This thus shows that the value of $\mathrm{uB}((i_1,j_1),(i_2,j_2))$ is completely determined by the vertices from $\Gamma^1 \cup \Gamma^2 \cup \cdots \cup \Gamma^{i_1-1} \cup \Gamma^{i_2+1} \cup \Gamma^{i_2+2} \cup \cdots \cup \Gamma^n$. This simple observation makes the computation of $\mathrm{uB}(\Gamma)$ much easier. To illustrate how this can be done for each $n$ and $m$ we provide some details for the cases $m = 3$ and $m = 4$. For ease of reference we state the following elementary fact which will be used several times in the corresponding proofs:
\begin{equation}
\label{eq:sum}
	\sum_{i=1}^{n-1} |n-2i| = \left\{\begin{array}{cc}
		\frac{n(n-2)}{2} & \text{if } n\ \text{is even},\\
		\frac{(n-1)^2}{2} & \text{if } n\ \text{is odd}.\end{array}\right.
\end{equation}

\begin{proposition}
\label{pro:PnC3}
Let $n$ be a positive integer. Then
$$
	\mathrm{uB}(P_n \square C_3) = \left\{\begin{array}{cc}
		\frac{3n(n-2)(6n-1)}{4} & \text{if } n\ \text{is even},\\
		\frac{3(n-1)(2n+1)(3n-5)}{4} & \text{if } n\ \text{is odd}.\end{array}\right.
$$
\end{proposition}

\begin{proof}
As mentioned in the discussion preceding the statement of the proposition the graphs $P_1 \square C_3 = C_3$ and $P_2 \square C_3$ are highly distance-balanced. Since the two expressions from the proposition for $n = 1$ and $n = 2$ indeed result in $0$, we can thus assume $n \geq 3$ for the rest of the proof. Denote $\Gamma = P_n \square C_3$.

Let us determine the value $\mathrm{uB}((i_1,j_1),(i_2,j_2))$ for each pair of vertices where $i_1 \leq i_2$, $j_1, j_2 \in \mathbb{Z}_3$. In view of the above mentioned automorphisms $\alpha$ and $\beta$ we only need to consider the possibilities with $i_1 < i_2$. Moreover, by the argument from the paragraph preceding this proposition we only need to consider the vertices $(i,j)$ with $i < i_1$ or $i > i_2$. Now, for any such vertex $(i,j)$ its distance to $(i_1,j_1)$ is equal to $|i-i_1|$ or $|i-i_1|+1$, depending on whether $j = j_1$ or not, respectively. Similarly, the distance from $(i,j)$ to $(i_2,j_2)$ is equal to $|i-i_2|$ or $|i-i_2|+1$, depending on whether $j = j_2$ or not, respectively. Therefore, unless $i_2 = i_1+1$ and $j_1 \neq j_2$, we see that all the columns $\Gamma^i$ with $i < i_1$ belong to $W_{(i_1,j_1),(i_2,j_2)}$ while all the columns $\Gamma^i$ with $i > i_2$ belong to $W_{(i_2,j_2),(i_1,j_1)}$. In this case we thus get
$$
	\mathrm{uB}((i_1,j_1),(i_2,j_2)) = 3|(n-i_2) - (i_1-1)| = 3|n+1-i_1-i_2| = 3\cdot \mathrm{uB}_{P_n}(i_1,i_2).
$$
If however $i_2 = i_1+1$ and $j_1 \neq j_2$, then each vertex of the form $(i,j_2)$ with $i < i_1$, and similarly each vertex of the form $(i,j_1)$ with $i > i_2$, has equal distance to $(i_1,j_1)$ and $(i_2,j_2)$. In this case we thus get
$$
	\mathrm{uB}((i_1,j_1),(i_1+1,j_2)) = 2|n-2i_1|.
$$
It thus follows that $\mathrm{uB}(\Gamma) = 3^3 \cdot \mathrm{uB}(P_n) - 3\cdot 2\cdot \sum_{i = 1}^{n-1}|n-2i|$, and so we can simply apply Proposition~\ref{pro:paths} and \eqref{eq:sum}.
\end{proof}

\begin{proposition}
\label{pro:PnC4}
Let $n$ be a positive integer. Then
$$
	\mathrm{uB}(P_n \square C_4) = \left\{\begin{array}{cc}
		\frac{2(n-2)(16n^2-13n+6)}{3} & \text{if } n\ \text{is even},\\
		\frac{2(n-1)(16n^2-29n+3)}{3} & \text{if } n\ \text{is odd}.\end{array}\right.
$$
\end{proposition}

\begin{proof}
The idea of the proof is similar as before. Denote $\Gamma =  P_n \square C_4$. We start with $4^3\cdot \mathrm{uB}(P_n)$ and then subtract an appropriate number to compensate for the pairs $(i_1,j_1), (i_2,j_2)$, $i_1 < i_2$, where it is not true that entire columns $\Gamma^i$ with $i < i_1$ are contained in $W_{(i_1,j_1),(i_2,j_2)}$. Of course, this time such a situation can occur in different ways. One is for sure if $i_2 = i_1+1$ and $j_1 \neq j_2$ (there are two essentially different possibilities here), but there is also the possibility that $i_2 = i_1 + 2$ and $j_2 = j_1 + 2$. Let us consider each of these three possibilities in some detail.

If $i_2 = i_1+1$ and $j_2 = j_1 \pm 1$, then for each $i < i_1$ the vertices of the form $(i,j)$ where $j \in \{j_2, j_1+2\}$ have equal distance to $(i_1,j_1)$ and $(i_2,j_2)$, while the remaining vertices of the form $(i,j)$ are in $W_{(i_1,j_1),(i_2,j_2)}$ (and similarly for the vertices $(i,j)$ with $i > i_2$). To compensate for these pairs we thus have to subtract $4\cdot 2 \cdot 2 \cdot \sum_{i=1}^{n-1}|n-2i|$ from $4^3\cdot \mathrm{uB}(P_n)$.

If $i_2 = i_1+1$ and $j_2 = j_1 + 2$, then for each $i < i_1$ the vertices of the form $(i,j_2)$ are in $W_{(i_2,j_2),(i_1,j_1)}$ while the remaining vertices of the form $(i,j)$ are in $W_{(i_1,j_1),(i_2,j_2)}$. To compensate for these pairs we thus have to subtract $4\cdot 1 \cdot 2 \cdot \sum_{i=1}^{n-1}|n-2i|$ from $4^3\cdot \mathrm{uB}(P_n)$.

Finally, if $i_2 = i_1 + 2$ and $j_2 = j_1 + 2$, then for each $i < i_1$ the vertices of the form $(i,j_2)$ have equal distance to $(i_1,j_1)$ and $(i_2,j_2)$ while the remaining vertices of the form $(i,j)$ are in $W_{(i_1,j_1),(i_2,j_2)}$. To compensate for these pairs we thus have to subtract $4\cdot 1 \cdot 1 \cdot \sum_{i=1}^{n-2}|(n-1)-2i|$ from $4^3\cdot \mathrm{uB}(P_n)$.

We can now apply Proposition~\ref{pro:paths} and \eqref{eq:sum} to obtain the stated result.
\end{proof}

Using a similar approach one can thus obtain the formula for $\mathrm{uB}(P_n \square C_m)$ for any particular choice of $m \geq 3$. We state just one more such result (without proof).

\begin{proposition}
\label{pro:PnC5}
Let $n$ be a positive integer. Then
$$
	\mathrm{uB}(P_n \square C_5) = \left\{\begin{array}{cc}
		\frac{5(n-2)(50n^2-47n+24)}{12} & \text{if } n\ \text{is even},\\
		\frac{5(n-1)(50n^2-97n+21)}{12} & \text{if } n\ \text{is odd}.\end{array}\right.
$$
\end{proposition}

\section{Distance-unbalancedness of trees}
\label{sec:trees}

In the previous section we determined the distance-unbalancedness index of members of a few well-known families of graphs. The stars $S_n$ and the paths $P_n$, which were among the studied families, are examples of trees. Observe that the distance-unbalancedness index of a star of order $n$ (which is $(n-1)(n-2)$) is much smaller than that of the path of the same order (which is roughly $n^3/6$). This is in stark contrast with the fact proved in~\cite{DMSTZ18} that the path of order $n$ has the smallest possible Mostar index among all trees of order $n$ while the star of order $n$ has the larges possible Mostar index among all trees of order $n$. This phenomenon can be explained to some extent by observing that the star of order $n$ has diameter $2$ and in fact its Mostar index coincides with its distance-unbalancedness, while on the other hand the path of order $n$ has diameter $n-1$, and its Mostar index forms just a tiny part of its distance-unbalancedness.

It is the aim of this section to initiate an investigation of distance-unbalancedness of trees. It seems that determining this index for all trees of order $n$ is a very difficult (if not impossible) task. We are thus primarily interested in obtaining some preliminary results, making some computer-assisted investigations and identifying some interesting problems that might be feasible to solve in the future.

Note that for any tree $\Gamma$ of order $n \geq 2$ and its leaf $u$ whose unique neighbor in $\Gamma$ is $v$ we have that $\mathrm{uB}(u,v) = n-2$. This shows that the only highly distance-balanced trees are the trivial ones $K_1$ and $K_2$. For all of the remaining ones the distance-unbalancedness index is nonzero. One of the most fundamental problems to consider is thus for sure the following one:

\begin{problem}
\label{pro:treemin}
For each integer $n \geq 3$ determine the minimum and the maximum value of $\mathrm{uB}(\Gamma)$ where $\Gamma$ is a tree of order $n$ and classify all trees that attain these two values.
\end{problem}

Using a computer one can of course determine the minimum and maximum value from the above problem at least for some small orders $n$. We investigated all trees up to order $15$ and obtained some interesting data. In the following table, for each order $n$, the smallest (min) and second smallest (min') possible values of distance-unbalancedness for trees of order $n$ are given, as well as the second largest (max') and largest (max) value of this index. We also state the number of all tress of given order ($\#\text{all}$), the number of trees attaining the smallest ($\#\text{min}$) and the number of trees attaining the largest ($\#\text{max}$) possible value of distance-unbalancedness.
$$
\begin{array}{c|ccccccccccccc}
	\text{n} & 		3 & 4 & 5   & 6   & 7  & 8 & 9 & 10 & 11 & 12 & 13 & 14 & 15\\
	\hline \hline
	\text{min} &  		2 & 6 & 12 & 20 &	30  & 42	& 56  &  72  &  90 &  110 & 132 & 156 & 182    \\
	\text{min'} &		  &    & 14 & 24 & 38  & 54 & 74  & 96  &  122 & 150 & 182 & 216 & 254 \\
	\text{max'} &         &    & 14 & 30 & 54  & 88 & 134 & 192 & 277 & 372 & 494 & 636 & 806 \\
	\text{max} &       2 & 6 & 16 & 32 &	56  & 90 & 138 &  198 & 278  & 378 & 495 & 638 & 808 \\ \hline
	\#\text{all} &		1 & 2 & 3   & 6  & 11  & 23	& 47  & 106 & 235  & 551 & 1301 & 3159 & 7741	\\
	\#\text{min} &     1 & 2 & 1   &  1 &	1    & 1  &   1  &  1   &  1 & 1  & 1 & 1 & 1 \\
	\#\text{max} & 1 & 2 & 1   & 	1 &	1   & 1   & 1    &  3  &  1 & 1 &  1 & 1 & 1 \\
\end{array}
$$
One can now analyze the trees attaining the smallest (few) and largest (few) values of distance-unbalancedness of given small orders. It turns out that for each $n$ with $5 \leq n \leq 15$ the unique tree $T$ of order $n$ with the smallest possible value of $\mathrm{uB}(T)$ among all trees of order $n$ is the star $S_n$ (with $\mathrm{uB}(S_n) = (n-1)(n-2)$). Moreover, it turns out that for each $n$ with $5 \leq n \leq 15$ there is a unique tree $T$ of order $n$ having the second smallest possible value of $\mathrm{uB}(T)$ among all trees of order $n$. This tree is isomorphic to the merged star $\SS{m}{m}$ with $n = 2m+2$ when $n$ is even, and to the subdivided merged star $\SSx{m}{m}$ with $n = 2m+3$ when $n$ is odd. Based on these observations we make the following conjectures.

\begin{conjecture}
\label{con1}
For any $n \geq 5$ the star $S_{n-1}$ is the unique tree of order $n$ having the smallest possible distance-unbalancedness index among all trees of order $n$.
\end{conjecture}

\begin{conjecture}
\label{con2}
For any $n \geq 5$, depending on whether $n$ is even or odd, respectively, the merged star $\SS{(n-2)/2}{(n-2)/2}$ or the subdivided merged star $\SSx{(n-3)/2}{(n-3)/2}$, respectively, is the unique tree of order $n$ having the second smallest possible distance-unbalancedness index among all trees of order $n$.
\end{conjecture}

It seems that proving the above two conjectures (especially the second one) will not be very easy. However, in the next few results we at least determine the distance-unbalancedness index of all merged stars and subdivided merged stars and show that for an even number $2n \geq 6$ the merged star $\SS{n-1}{n-1}$ is the unique tree among all merged stars and subdivided merged stars of order $2n$ that has the smallest possible distance-unbalancedness index. Similarly, we show that for an odd number $2n+1 \geq 5$ the subdivided merged star $\SSx{n-1}{n-1}$ is the unique tree among all merged stars and subdivided merged stars of order $2n$ that has the smallest possible distance-unbalancedness index.

\begin{proposition}
	\label{dstar:prop1}
	Let $n$ and $m$ be positive integers with $n \ge m \ge 1$. Then
	$$
	  \mathrm{uB}(\SS{n}{m})=(n+m)^2+(n-m)(nm+1)+2nm.
	$$
\end{proposition}
\begin{proof}
	Denote $\Gamma = \SS{n}{m}$ and observe that since $n \ge m \ge 1$, the diameter of $\Gamma$ is 3. It is now easy to see that
	$$
	  \mathrm{Mo}^1(\Gamma) = (n+m)^2+(n-m), \quad
	\mathrm{Mo}^2(\Gamma) = 2nm\quad \text{and}\quad
	\mathrm{Mo}^3(\Gamma) = nm(n-m),
	$$
and so the result follows.
\end{proof}

\begin{proposition}
	\label{dstar:prop2}
	Let $n,m$ be positive integers with $n \ge m \ge 1$. Then
	$$
	\mathrm{uB}(\SSx{n}{m})=(n+m)(n+m+1)+(n+1)|n-m-1|+(m+1)(3n-m+1)+nm(n-m).
	$$
\end{proposition}
\begin{proof}
Observe first that the diameter of $\Gamma = \SSx{n}{m}$ is 4. It is easy to see that
	$$
	\mathrm{Mo}^1(\Gamma) = (n+m)(n+m+1)+|n-(m+1)| + n-m+1,
	$$
	$$
	\mathrm{Mo}^2(\Gamma) = n(m+1) + m(n+1) + n-m,
	$$
$$
	\mathrm{Mo}^3(\Gamma) = n|n-(m+1)| + m(n+1-m) \quad \text{and}\quad \mathrm{Mo}^4(\Gamma) = nm(n-m),
$$
and so the result follows.
\end{proof}

\begin{theorem}
\label{thm:mstars}
Let $N \geq 3$ be an integer and let $\mathcal{S}_N$ be the set of all merged stars and all subdivided merged stars of order $N$. Then the following hold:
\begin{itemize}
\itemsep = 0pt
\item[(i)] If $N = 2(m+1)$ is even, then $\SS{m}{m}$ is the unique member of $\mathcal{S}_N$ having the smallest possible distance-unbalancedness index among all members of $\mathcal{S}_N$ and $\mathrm{uB}(\SS{m}{m}) = 6m^2$.
\item[(ii)] If $N = 2m+3$ is odd, then $\SSx{m}{m}$ is the unique member of $\mathcal{S}_N$ having the smallest possible distance-unbalancedness index among all members of $\mathcal{S}_N$ and $\mathrm{uB}(\SSx{m}{m}) = 6m^2+6m+2$.
\end{itemize}
\end{theorem}

\begin{proof}That the result holds for all $N$ with $3 \leq N \leq 15$ can easily be verified using Propositions~\ref{dstar:prop1} and~\ref{dstar:prop2}. For the rest of the proof we thus assume $N \geq 16$.

Consider first all merged stars $\SS{n}{m}$ where $N = n+m+2$ with $n \geq m \geq 1$. Of course, $m = N-n-2$ and $1 \leq m \leq (N-2)/2$ while $(N-2)/2 \leq n \leq N-3$. Write $f(n) = \mathrm{uB}(\SS{n}{m})$. Using Proposition~\ref{dstar:prop1} and rearranging terms we find that
$$
	f(n) = -2n^3 + (3N-8)n^2 -(N^2-6N+6)n + (N-2)(N-3)
$$
is a polynomial in $n$ of degree $3$ with a negative leading coefficient. As $f(0) > 0$ and $f(2) = -N^2+19N-54 < 0$ (recall that $N \geq 16$), $f$ has a root in the interval $(0,2)$. Since $f((N-8)/2) = 30$ and $(N-8)/2 \geq 4 > 2$, $f$ has another root in the interval $(2,(N-8)/2)$. Finally, $f(N-2)=N^2-3N+2 > 0$ and $f(N-1)=6-4N < 0$, and so $f$ has its third root in the interval $(N-2,N-1)$. Since we are only interested in the value of $f(n)$ for $n$ satisfying $(N-2)/2 \leq n \leq N-3$, we thus only need to compare $f(N-3)$ to $f(n)$ for the smallest integer $n$ with $n \geq (N-2)/2$. Now,
$$
	f(N-3) = 2(N-1)(N-3), \ f((N-2)/2) = \frac{3(N-2)^2}{2}\ \text{and}\ f((N-1)/2) = \frac{7N^2-28N+29}{4},
$$
and so $f(N-3) > f((N-1)/2) > f((N-2)/2)$, implying that for $N$ even $\mathrm{uB}(\SS{n}{m})$ is minimal for $m = n = (N-2)/2$ (and equals $3(N-2)^2/2 = 6n^2$), while for $N$ odd $\mathrm{uB}(\SS{n}{m})$ is minimal for $n = (N-1)/2$ and $m = (N-3)/2$ (and equals $(7N^2-28N+29)/4$).

Consider now all subdivided merged stars $\SSx{m}{n}$ where $N = n+m+3$ with $n \geq m \geq 1$. Then $m = N-n-3$, $1 \leq m \leq (N-3)/2$ and $(N-3)/2 \leq n \leq N-4$. Write $g(n) = \mathrm{uB}(\SSx{n}{m})$. If $m = n$ (which is only possible if $N = 2n+3$ is odd), then Proposition~\ref{dstar:prop2} implies that
$$
	g(n) = 6n^2+6n+2 = \frac{3N^2-12N+13}{2}.
$$
If however $n \geq m+1$ (in which case $n \geq (N-2)/2$), then Proposition~\ref{dstar:prop2} implies that
$$
	g(n) = -2n^3 + (3N-11)n^2 - (N^2-10N+17)n
$$
is a polynomial in $n$ of degree $3$ with a negative leading coefficient. Note that $g(0) = 0$. Similarly as above we find that $g$ also has one root in the interval $((N-10)/2,(N-8)/2)$ and another one in the interval $(N-2,N-1)$. Therefore, in the interval $((N-2)/2, N-4)$ the function $g$ attains its minimum either in $N-4$ or in the smallest integer $n$ with $n \geq (N-2)/2$. As
$$
	g(N-4) = (N-4)(3N-5),\ g((N-2)/2) = \frac{7N^2-30N+32}{4}\ \text{and}\ g((N-1)/2) = 2(N-1)(N-3),
$$
we see that $g(N-4) > g((N-1)/2) > g((N-2)/2)$ (recall that $N \geq 16$), and so the fact that also $(3N^2-12N+13)/2 < 2(N-1)(N-3)$ implies that for $N$ odd $\mathrm{uB}(\SSx{n}{m})$ is minimal for $m = n = (N-3)/2$ (and equals $(3N^2-12N+13)/2 = 6n^2+6n+2$), while for $N$ even $\mathrm{uB}(\SS{n}{m})$ is minimal for $n = (N-2)/2$ and $m = (N-4)/2$ (and equals $(7N^2-30N+32)/4$).

To complete the proof we now simply have to compare the above obtained smallest possible value of $\mathrm{uB}(\SS{n}{m})$ and $\mathrm{uB}(\SSx{n}{m})$ depending on whether $N$ is even or odd.
\end{proof}

Returning to the above data regarding the distance-unbalancedness index of trees up to order $15$, it seems that determining the largest (or second largest) possible value of this index among all trees of order $n$ is even more difficult than proving Conjecture~\ref{con1}. Inspecting the trees attaining these values for $n$ up to $15$ one finds that all of these trees are spider graphs. As was indicated in the above table, for each integer $n$ with $5 \leq n \leq 15$, except for $n = 10$, there is a unique tree of order $n$ having the largest possible distance-unbalancedness index among all trees of order $n$. For $n = 10$ there are three such trees, namely the spiders $\mathrm{Sp}(3,2,2,2), \mathrm{Sp}(3,3,2,1)$ and $\mathrm{Sp}(3,2,2,1,1)$. For $n$ with $5 \leq n \leq 15$, $n \neq 10$, the unique maximal tree with respect to the distance-unbalancedness index is given in the following table, where for $\mathrm{Sp}(n_1, n_2, \ldots , n_k)$ we simply state the sequence $(n_1,n_2, \ldots , n_k)$ (see also Figure~\ref{fig:spiders}).

{\small
\begin{tabular}{c|c@{\,}c@{\,}c@{\,}c@{\,}c@{\,}c@{\,}c@{\,}c@{\,}c@{\,}c}
	n & 		5   & 6   & 7  & 8 & 9 & 11 & 12 & 13 & 14 & 15 \\
	\hline
	max &  (2,1,1) & (2,2,1) & (2,2,1,1) & (2,2,2,1) & (3,2,2,1) & (3,3,2,2) & (3,3,2,2,1) & (3,3,3,2,1) & (3,3,3,2,2) & (4,3,3,2,2)
\end{tabular}
}
\begin{figure}[!h]
\begin{center}
	\includegraphics[scale=0.9]{max_spiders.pdf}
	\caption{Maximal trees with respect to the distance-unbalancedness index.}
	\label{fig:spiders}
\end{center}
\end{figure}
\medskip

In view of this data it is interesting to ask whether for all $n \geq 11$ there is a unique tree with maximal distance-unbalancedness index among all trees of order $n$ or there are some values of $n$ where a similar situation as for $n = 10$ occurs. Also, the following problem is of interest.

\begin{problem}
\label{prob:spiders}
For each integer $n \geq 5$ classify the spiders of order $n$ that have the largest possible distance-unbalancedness index among all spiders of order $n$.
\end{problem}

\section{Minimal distance-unbalancedness}
\label{sec:min}

In the previous section we considered the problem of determining the smallest nonzero and the largest possible distance-unbalancedness index among all trees of given order. Given that those two problems appear to be rather difficult it seems that determining the smallest nonzero (and the largest) possible distance-unbalancedness index among all connected graphs of given order might be even more difficult. It seems natural to speculate that, as with trees, the smallest possible nonzero distance-unbalancedness index should increase with the order of the considered graphs. Nevertheless, as we shall see later, this is not true.

One can of course again use computers to investigate all connected graphs up to a given reasonably small order. We have looked at all connected graphs up to order $10$ using the lists of all such graphs available at~\cite{McKayWeb}. For each order $n$, $3 \leq n \leq 10$, we determined the smallest possible nonzero value $\mathrm{uB}(\Gamma)$ among all connected graphs of order $n$ and determined all graphs $\Gamma$ of order $n$ attaining this value. The following table summarizes the obtained data where for each order $n$ the smallest possible nonzero distance-unbalancedness index (min) for connected graphs of order $n$ and the number of graphs attaining this value ($\#\text{min}$) are given.
$$
\begin{array}{c|cccccccc}
	\text{n} & 		3 & 4 & 5   & 6   & 7  & 8 & 9 & 10 \\
	\hline
	\text{min} &  		2 & 4 & 4 & 8   &	6  & 8  & 6  & 8 \\
	\#\text{min} &     1 & 1 & 2 & 10  &	8  & 2  & 1  & 17 \\
\end{array}
$$

With trees the computer generated data is large enough to give a very strong indication of what should be proved (Conjecture~\ref{con1}). In this general setting the set of orders considered is perhaps a bit to small to draw conclusions but it does indicate that perhaps the above mentioned ``natural speculation'' is not correct. Nevertheless, we propose the following problem (although the obtained data suggests that at least the second part of the problem might be extremely difficult if not impossible to solve).

\begin{problem}
\label{prob:allmin}
For each integer $n \geq 3$ determine the smallest possible nonzero value of the distance-unbalancedness index among all connected graphs of order $n$ and classify all graphs attaining this value.
\end{problem}

Looking at the values from the above table it seems that perhaps the solution to the first part of the problem is that for $n \geq 6$ the minimal nonzero value of the distance-unbalancedness index among all connected graphs of order $n$ is $6$ or $8$, depending on whether $n$ is odd or even, respectively. Of course, the obtained data is much too small to justify this speculation. Trying to determine whether it is indeed plausible that there does exist a constant $C$ such that for each integer $n$ the smallest nonzero distance-unbalancedness index among all connected graphs of order $n$ is at most $C$, it makes sense to have a closer look at the minimal graphs obtained for orders up to $n = 10$ (those of order $n$ that have the minimal possible nonzero distance-unbalancedness index among all connected graphs of order $n$).

Doing so we discovered two graphs that fall into a known infinite family. Namely, the unique minimal graph of order $9$ and one of the 10 minimal graphs of order $6$ belong to an infinite family of bi-regular graphs (some vertices are of valence $3$ while others are of valence $4$) constructed in~\cite{IlKlMi10} (see Figure~\ref{fig:kites} where these two graphs and the members of the family of orders $12$ and $15$ are illustrated). This suggests that determining the distance-unbalancedness index of all members of this infinite family might be of interest. Interestingly enough, the graphs from the family are of two very different types regarding $\ell$-distance-balancedness. Those of even order are not $1$-distance-balanced but are $\ell$-distance-balanced for all $\ell \geq 2$, while those of odd order are $1$-distance-balanced but are not $\ell$-distance-balanced for any odd $\ell \geq 3$.
\begin{figure}[!h]
\begin{center}
	\includegraphics[scale=0.75]{kite_graphs.pdf}
	\caption{The kite graphs $\mathrm{Ki}(n)$ for $2 \leq n \leq 5$.}
	\label{fig:kites}
\end{center}
\end{figure}

\begin{construction}
For each integer $n \geq 2$ the {\em Kite graph} $\mathrm{Ki}(n)$ is the graph of order $3n$ with vertex set
$$
	\{a_i \colon i \in \mathbb{Z}_n\} \cup \{b_i \colon i \in \mathbb{Z}_n\} \cup \{c_i \colon i \in \mathbb{Z}_n\},
$$
where for each $i \in \mathbb{Z}_n$ the vertices $b_i$ and $c_i$ are adjacent to each other and to $a_i$ and $a_{i+1}$ (the indices are computed modulo $n$). The reason we chose this name for them is that we can think of obtaining $\mathrm{Ki}(n)$ by taking $n$ copies of the $4$-cycle with one additional chord (which can be viewed as a kite) and then gluing them along vertices of valence $2$ in a cyclic fashion (so that the ``old'' vertices of valence $2$ become of valence $4$ in the obtained graph).
\end{construction}

\begin{proposition}
\label{pro:zmaji}
Let $n \geq 2$ be an integer. Then the graph $\mathrm{Ki}(n)$ is a bi-regular graph of diameter $n$ with $n$ vertices of valence $4$ and $2n$ vertices of valence $3$. Moreover, the following hold:
\begin{itemize}
\itemsep = 0pt
\item[(i)] If $n$ is even, then $\mathrm{Ki}(n)$ is not $1$-distance-balanced but is $\ell$-distance-balanced for all $\ell$, $2 \leq \ell \leq n$. Moreover, $\mathrm{uB}(\mathrm{Ki}(n)) = 4n$.
\item[(ii)] If $n$ is odd, then $\mathrm{Ki}(n)$ is $\ell$-distance-balanced for $\ell=1$ and all even $\ell$ with $\ell \leq n-1$, but is not $\ell$-distance-balanced for any odd $\ell$, $3 \leq \ell \leq n$. Moreover, $\mathrm{uB}(\mathrm{Ki}(n)) = 2n(n-2)$.
\end{itemize}
\end{proposition}

\begin{proof}
The initial assertions follow straight from the definition. Let $\rho$ and $\sigma$ be the permutations of the vertex-set of $\mathrm{Ki}(n)$ such that for each $i \in \mathbb{Z}_n$
$$
	a_i\rho = a_{i+1}, b_i\rho = b_{i+1}, c_i\rho = c_{i+1},\ \text{and}\ a_i\sigma = a_{-i}, b_i\sigma = b_{-(i+1)}, c_i\sigma = c_{-(i+1)}.
$$
Furthermore, le $\tau_0$ be the permutation of the vertex-set of $\mathrm{Ki}(n)$, interchanging $b_0$ with $c_0$ and leaving all other vertices fixed. It is easy to see that these three permutations are in fact automorphisms of $\mathrm{Ki}(n)$. The automorphism group of $\mathrm{Ki}(n)$ thus contains the group $\langle \rho, \sigma, \tau_0\rangle$, and so it is clear that the only possible unbalanced pairs are of the form $\{a_i, b_j\}$ and $\{a_i, c_j\}$. Moreover, it suffices to determine all $\mathrm{uB}(a_0, b_j)$ where $j < n/2$. To do this observe that $c_0$ is the unique vertex which has equal distance to $a_0$ and $b_0$ while for $1 \leq j < n/2$ no vertex of $\mathrm{Ki}(n)$ has equal distance to $a_0$ and $b_j$. It is now straightforward to verify that if $n$ is even $\mathrm{uB}(a_0,b_0) = 1$ and $\mathrm{uB}(a_0,b_j) = 0$ for all $j$, $1 \leq j < n/2$, while if $n$ is odd $\mathrm{uB}(a_0,b_0) = 0$ (see also~\cite{IlKlMi10}) and $\mathrm{uB}(a_0, b_j) = 1$ for all $j$, $1 \leq j < n/2$. To complete the proof one now simply has to count the number of all such pairs. We leave details to the reader.
\end{proof}

The above proposition shows that the graph $\mathrm{Ki}(n)$, where $n$ is odd, has a rather large distance-unbalancedness index, and so it is not surprising that there exist graphs of order $3n$ with a smaller nonzero distance-unbalancedness index (except of course for $n = 3$). Indeed, for any odd $n \geq 5$, Proposition~\ref{pro:multi} and Proposition~\ref{pro:zmaji} show that $\mathrm{uB}(\mathrm{Ki}(n)) = 2n(n-2) \geq 6n > 6n-4 = \mathrm{uB}(\Gamma)$, where $\Gamma$ is the complete multipartite graph of order $3n$ of the form $K_{2,1,1,\ldots , 1}$.

While on the other hand the graphs $\mathrm{Ki}(n)$ with $n$ even seem to be better in this respect their distance-unbalancedness index is still not bounded by a constant. But there is a different generalization of the graph $\mathrm{Ki}(2)$ (which is minimal among connected graphs of order $6$ in the above mentioned sense) that gives such a constant at least for graphs of even orders. Observe that (but see also~\cite{IlKlMi10}) the construction of the graph $\mathrm{Ki}(n)$ can also be described as follows. Start with the cycle $C_{2n}$ of length $2n$ and think of its vertices being labeled by elements of $\mathbb{Z}_{2n}$ in a natural way. For each even $i \in \mathbb{Z}_{2n}$ introduce a copy $i'$ of the vertex $i$ and connect it to each of $i$, $i-1$ and $i+1$. The resulting graph is of course $\mathrm{Ki}(n)$. But what if we introduce copies of just a few vertices and not of every other one in $C_{2n}$?

\begin{construction}
Let $n \geq 2$ be an integer. Label the vertices of the cycle $C_{2n}$ in a natural way by elements of $\mathbb{Z}_{2n}$ so that for each $i$ the vertex $i$ is adjacent to $i+1$ and $i-1$. Introduce two new vertices, say $x$ and $y$, and connect $x$ to $0$, $-1$ and $1$, and $y$ to $n$, $n-1$ and $n+1$. Denote the obtained graph of order $2n+2$ by $\tilde{C}_{2n}$ (see Figure~\ref{fig:cycles2} where the first few examples are illustrated).
\end{construction}
\begin{figure}[!h]
\begin{center}
	\includegraphics[scale=0.75]{cycles2.pdf}
	\caption{The graphs $\tilde{C}_{2n}$ for $2 \leq n \leq 5$.}
	\label{fig:cycles2}
\end{center}
\end{figure}

\begin{proposition}
\label{pro:CQ}
Let $n \geq 2$ be an integer. Then $\mathrm{uB}(\tilde{C}_{2n}) = 8$. Moreover, the graph $\tilde{C}_{2n}$ is of order $2n+2$ and diameter $n$ and is not $1$-distance-balanced, but is $\ell$-distance-balanced for all $\ell$, $2 \leq \ell \leq n$.
\end{proposition}

\begin{proof}
Let $\tau_x, \tau_y$ be the automorphisms of $\tilde{C}_{2n}$ fixing all of its vertices except $0$ and $x$, respectively $n$ and $y$, which are swapped. Similarly, let $\tau_0$ be the automorphism interchanging vertices $i$ and $-i$ of $\tilde{C}_{2n}$ for all $i \in \mathbb{Z}_{2n}$ and fixing each of $0$, $n$, $x$ and $y$. Moreover, let $\rho$ be the automorphism interchanging each $i \in \mathbb{Z}_{2n}$ with $n+i$ and interchanging $x$ with $y$. In view of these automorphisms it thus suffices to consider only pairs of vertices of the form $\{i,j\}$, where $0 \leq i < j \leq n-1$. In fact, because of the automorphism $\tau_0\rho$ we can in fact assume that $i < n-j$ (if $i = n-j$ this automorphism actually swaps $i$ and $j$, and so this pair is balanced). Now, if $i \neq 0$ or $j \neq 1$, then $x \in W_{i,j}$ while $y \in W_{j,i}$, and so the pair $\{i,j\}$ is clearly balanced. On the other hand, as $x$ is a common neighbor of $0$ and $1$, the pair $\{0,1\}$ is not balanced and $\mathrm{uB}(0,1) = 1$ (as $y \in W_{1,0}$). Therefore, the only unbalanced pairs correspond to the 8 edges $\{0, \pm 1\}$, $\{x, \pm 1\}$, $\{n, n \pm 1\}$, $\{y, n \pm 1\}$ and so the result follows.
\end{proof}

Perhaps there is a similar construction (based on a suitable cycle) that for each odd $n \geq 7$ results in a graph of order $n$ with distance-unbalancedness index equal to $6$, but it is not completely obvious how to obtain it. We thus leave the question of whether (and how) such graphs might be constructed unanswered.
\bigskip

We conclude the paper with the following remarks. Even though we were able to obtain a connected graph $\Gamma$ of even order $n$ with $\mathrm{uB}(\Gamma) = 8$ for each even $n \geq 6$, it does feel that in obtaining these graphs we were ``cheating'' a bit - the examples were obtained by taking a cycle (which is in fact distance transitive, implying that for any pair of vertices there is an automorphism swapping them) and then changing it very slightly so as to introduce a tiny unbalancedness. But what if one only considers regular graphs? It is well known and easy to see that all (connected) regular graphs up to order $6$ are Cayley graphs of abelian groups and are thus highly distance-balanced. The only connected regular graph of order $7$ which is not a Cayley graph (of an abelian group) is the complement of the disjoint union of a $3$-cycle and a $4$-cycle. It is not difficult to see however, that this graph is highly distance-balanced.

With order $8$ the situation changes. In fact, the two connected graphs of order $8$ with minimal nonzero distance-unbalancedness index (recall that it is equal to $8$) are both cubic. One of them is $\tilde{C}_6$, while the other is obtained from the $8$-cycle by adding the edges $\{0,4\}, \{1,3\}, \{2,6\}$ and $\{5,7\}$. Using a computer one can verify that all connected regular graphs of order $9$ are again highly distance-balanced. For order $10$ recall that the minimal nonzero distance-unbalancedness index among all connected graphs of order $10$ is $8$ and that it is attained by $17$ graphs. It turns out that apart from $\tilde{C}_8$ (which is not regular) all of the remaining $16$ graphs are $4$-regular. Finally, using a computer one can verify that the smallest nonzero distance-unbalancedness index among all connected regular graphs of order $11$ is $10$ and that it is attained by $43$ graphs.

It is not difficult to see that $9$ is the largest odd order $n$ such that there does not exist a connected regular graph of order $n$ which is not highly distance-balanced. Namely, let $n = 2m+1$ with $m \geq 5$. Take any connected $4$-valent graph $\Delta$ of order $m$ and let $x$ and $y$ be a pair of adjacent vertices in $\Delta$. Construct a new graph $\Gamma$ by taking two copies of $\Delta$ (denote the second copy by $\Delta'$ and the analogues of $x$ and $y$ in $\Delta'$ by $x'$ and $y'$), deleting the edges $xy$ and $x'y'$, introducing a new vertex $z$ and connecting it to $x, y, x'$ and $y'$. The obtained graph is clearly a connected $4$-regular graph of order $n$, but the pair $\{z,x\}$ is not balanced as $W_{z,x}$ contains at least $z$, $y$ and the whole copy $\Delta'$. In a similar way a connected cubic graph of order $2n$ for each $n \geq 5$ which is not highly distance-balanced can be constructed (this time start with two copies of a cubic graph of order $n-1$ and then add a pair of connected vertices in a similar way as before, but see also~\cite[Theorem~5.7]{MikSpa18}). The following problem thus might be of interest.

\begin{problem}
\label{prob:reg}
For each positive integer $n \geq 10$ determine the smallest possible nonzero value of the distance-unbalancedness index among all connected regular graphs of order $n$ and classify all graphs attaining this value.
\end{problem}

\section*{Acknowledgements}

\noindent
\v S. Miklavi\v c acknowledges support by the Slovenian Research Agency (research program P1-0285 and research projects N1-0062, J1-9110, J1-1695, N1-0140, N1-0159, J1-2451).

\noindent
P.~\v Sparl acknowledges support by the Slovenian Research Agency (research program P1-0285 and research projects J1-9108, J1-9110, J1-1694, J1-1695, J1-2451).


\section{Introduction}
\label{intro}
The Boltzmann (transport) equation is an integro-differential equation which describes the local density of particles traveling and interacting inside a media~\cite{Duderstadt,CaseZweifel1967}. In high energy density physics phenomena, it appears in radiative transfer and describes the absorption and emission of photons within physical materials~\cite{pomraning2005equations}. Since the radiation opacity and the heat capacity of materials depend strongly on the material's temperature, it generates nonlinear waves, i.e. Marshak waves. The physics that is governed by the radiative transfer equation (RTE) is used for the study of many physical phenomena, for example, modeling theoretical astrophysics (supernova)~\cite{castor2004,pomraning2005equations}, laboratory plasma and inertial confinement fusion (ICF)~\cite{lindl2004,rosen1996}.

There are several ways to solve the transport equation, among them are the spherical harmonics method (the $P_N$ approximation)~\cite{pomraning2005equations}, the discrete ordinates method (the $S_N$ method)~\cite{sn} and Monte-Carlo simulations~\cite{imc}. Of course, the Monte-Carlo method is a statistical model while the first two solvers are deterministic transport equations and tend to the exact solution only for infinite set of coupled equations ($N\to\infty$). However, in some cases solving exactly the transport equation is impractical and resource-consuming. Thus, different approximations were introduced to reduce the complexity of the problem. The most well-known approximation, the {\em classic diffusion approximation}, can be developed by assuming that the angular dependence is isotropic or close to isotropic, and it is a simplification of the $P_N$ using $N=1$, (assuming steady-state behavior of the particle current)~\cite{pomraning2005equations}. Thus, the diffusion approximation tends to the exact solution in the limit of local thermodynamic equilibrium (LTE), close to the steady-state. As a result, one pitfall of the diffusion approximation is that it fails to describe the particles' density in highly anisotropic scenarios, specifically in optically thin media~\cite{su1997analytical,Olson1999}. Flux-Limiters (in a diffusion-like equation) or Variable Eddington Factors (in full $P_1$ equations) may be introduced, however these approximations are nonlinear making the solution hard to obtain for multi-dimensional cases, and yield unstable solutions in specific situations~\cite{Su2001}.

In the $P_N$ approximation the specific intensity is expanded in a full set of spherical harmonics thus, transforming the transport equation to a set of coupled equations, with a cutoff method or closure~\cite{Duderstadt,CaseZweifel1967,pomraning2005equations}. For example, the classical $P_N$ approximation is obtained by introduction of a closure that sets to zero all the coefficients of the expansion with an index larger than $N$~\cite{Duderstadt,CaseZweifel1967,pomraning2005equations}. Alternative $P_N$ closure was introduced by Pomraning~\cite{pomraning1964generalized, pomraning1965asymptotically} for time-independent problems. It was derived as a general-$N$ expansion of the asymptotic diffusion approximation~\cite{case1953introduction}, based on the spatial asymptotic time-independent solution for the exact Boltzmann equation. The physical justification for the asymptotic distribution of particles to fit well in various finite systems with sufficient number of mean free paths (MFP) comes from the physics of neutron transport. The asymptotic diffusion yields the correct critical radii, even for systems with less than two mean MFPs, while the classical $P_N$ approximation converges slowly compared to the exact solution~\cite{case1953introduction, BellGlasstone}.

For time-dependent problems, a related approximation was proposed by Heizler~\cite{heizler2010asymptotic}, namely the {\em asymptotic $P_1$ approximation} that is based on an asymptotic behavior in both space and time of the transport equation. This approximation is, in some sense, a time-dependent equivalent to the asymptotic diffusion approximation. This approximation has a $P_1$ form that overcomes the parabolic nature of the asymptotic diffusion. It yields the correct time-independent eigenvalue of the exact transport equation, and produces velocity that is very close to the correct one. In homogeneous problems, the asymptotic $P_1$ yields better results than the asymptotic diffusion, and much better than the classic $P_1$ approximation. This approximation was derived and tested also for local LTE radiative transfer~\cite{heizler2012}, and shares similar asymptotic behavior as $P_2$ approximation~\cite{heizler2012sp}. Further derivation for heterogeneous problems was introduced by the discontinuous asymptotic $P_1$ approximation~\cite{cohen2018discontinuous, cohen2019discontinuous}, which is important especially in radiative transfer problems, when the cross-sections are strong functions of the material's temperature.

Recently a new time-dependent asymptotic $P_N$ approximation was introduced~\cite{harel2020time}, which is a generalization of the work with the asymptotic $P_1$, and is based on the time-independent closure proposed by Pomraning while using time-asymptotic analysis. This approximation was tested in homogeneous media problems and yields very good results, including the {\em{tails}} of the distribution, even for low-order $N$. In this paper we generalize the testing to heterogeneous problems, specifically for radiative transfer and utilizing the Su-Olson benchmark~\cite{su1997analytical}. To avoid redundancy with~\cite{harel2020time}, we present only short summary of the derivation of the time-dependent asymptotic $P_N$ approximation, with  attention to the special modifications that were done to include radiative transfer problems.

First, in Section~\ref{section_pn} we present the $P_N$ for  both the classic and the Pomraning closure methods. Afterwards we derive our approximation named the time-dependent asymptotic $P_N$ approximation in its radiative transfer form. Then, in Sec.~\ref{test_benchmark} we compare our method to the classical $P_N$ approximation and the Su-Olson benchmark~\cite{su1997analytical}.

\section{$P_N$ Equations} \label{section_pn}

In this section, we present a derivation of the time-dependent asymptotic $P_N$ approximation in its radiative transfer form. First, we present the simple form of the $P_N$ equations with the classic $P_N$ closure, and the time-independent asymptotic closure proposed by Pomraning~\cite{pomraning1964generalized}. Next, a new time-dependent $P_N$ approximation, that is based on the a time-asymptotic analogy of Pomraning's closure, is described. This derivation is based on~\cite{harel2020time}, while taking special care of the new modifications that have to be made for radiative transfer.

\subsection{The Classic $P_N$ Approximation}

In one-dimensional slab geometry, the gray radiative transfer Boltzmann equation attains the form:
\begin{equation} \label{transport_Eq}
    \frac{1}{c}\pdv{I(x, \mu, t)}{t} + \mu \pdv{I(x, \mu, t)}{x} + \sigma_t I (x, \mu, t) = \sigma_a B(T_m(x,t))+ \frac{1}{2}\int_{-1}^{1} \sigma_s(\mu_0)I(x, \mu', t) d\mu' + S(x, t)
\end{equation}
Where $I(x, \mu, t)$ is the specific intensity of radiation in position $x$ at time $t$ in $\mu$ direction, i.e., the cosine of the angle with respect to the $x$-axis ($\mu=\Omega \cdot \hat{x}$). $c$ is the speed of light, $B(T_m(x,t))$ is the local black-body radiation emitted from the material, where $T_m(x,t)$ is the material temperature. $\sigma_t=\sigma_s+\sigma_a$ is the total cross section which is composed of $\sigma_s(\mu_0)$ (where $\mu_0\equiv\Omega \cdot \Omega'$), the scattering cross section (such as Thomson or Compton scattering) and $\sigma_a$, the absorbing cross section (the opacity). $S(x, t)$ is the external source term. The equation for the radiation is coupled to the complementary equation for the material temperature, $T_m(x,t)$:
\begin{equation}\label{materialequation}
    \frac{C_v(T_m(x,t))}{c} \pdv{T_m(x,t)}{t} = \sigma_a \left(\frac{1}{c}\int_{-1}^{1} I(x, \mu, t)d\mu - a T_m^4(x,t)\right)
\end{equation}
where $C_v(T_m(x,t))$ is the material's heat capacity and $a$ is the radiation constant and defined by $a\equiv 4\sigma_{\mathrm{sb}}/c$, and $\sigma_{\mathrm{sb}}$ is the Stefan-Boltzmann constant ($B(x,t) = aT_m^4$).

To derive the various $P_N$ approximations, we expand the specific intensity ($I(x,\mu,t)$) in a full set of spherical harmonics~\cite{Duderstadt,CaseZweifel1967,pomraning2005equations}, or in the one-dimensional case to a full set of Legendre polynomials:
\begin{subequations}
    \begin{equation}\label{polynomseries}
        I(x, \mu, t) = \sum_{n=0}^{\infty}\frac{2n+1}{2} I_n(x, t) P_n(\mu),
    \end{equation}
    \begin{equation}
    I_n(x,t) = \int_{-1}^{1} P_n(\mu) I(x, \mu, t)d\mu
\end{equation}
\end{subequations}
where $P_n(\mu)$ are the Legendre polynomials and the angular moments are given by $I_n(x,t)$. The first two moments are $I_0(x, t)\equiv cE(x, t)$ - the energy density and $I_1(x, t)\equiv F(x, t)$ - the radiation flux. Plugging Eq.~\ref{polynomseries} back into the transport equation (Eq.~\ref{transport_Eq}) and using the orthogonal property of the Legendre polynomials yields infinite coupled set of equations where each equation is given by:
\begin{equation}\label{pnequations}
\frac{1}{c} \pdv{I_n(x, t)}{t} + (\sigma_t-\sigma_s^{(n)}) I_n(x, t) +\left(\frac{n+1}{2n+1}\right)  \pdv{I_{n+1}(x, t)}{x} + \left(\frac{n}{2n+1}\right) \pdv{I_{n-1}(x, t)}{x} = S^{(n)}(x, t)
\end{equation}
$\sigma_s^{(n)}$ is the $n$'th moment of the scattering cross-section and $S^{(n)}(x, t)$ is the $n$'th moment of the external source term. Both of these terms can be derived in the same way that $I_n(x,t)$ is derived, i.e by expanding them into Legendre series.

In order to reduce the infinite sum to a finite sum some sort of closure needs to be introduced. As an example, for the classical $P_N$ approximation the closure is simply setting all the coefficients with index larger than $N$ to zero:
\begin{equation}
    \pdv{I_{n+1}(x, t)}{x}=0 \qquad\qquad n\geqslant N.
\end{equation}
This results in $N-1$ \textbf{exact} equations, that can be expressed with Eqs.~\ref{pnequations} and the last $N$'th \textbf{approximated} equation is defined by:
\begin{equation}\label{pn_classic_closure}
    \frac{1}{c}\pdv{I_N(x,t)}{t} + (\sigma_t - \sigma_s^{(N)})I_N(x,t) + \left(\frac{N}{2N+1}\right)\pdv{I_{N-1}(x, t)}{x} = S^{(N)}(x,t).
\end{equation}
Eqs.~\ref{pn_classic_closure} and~\ref{pnequations} yield a closed set of $N+1$ equations for $N+1$ unknown moments. The case of $N=1$ tends to the gray diffusion approximation, when neglecting the derivative of $F(x, t)$ (the radiation flux) with respect to time. The derivation, from here on, will use the notation of isotropic scattering term. For the case of anisotropic scattering, one have to replace $\sigma_t$ with $\sigma_t - \sigma_s^{(n)}$.

\subsection{Pomraning's Time-independent Asymptotic $P_N$}\label{pomramingpna}
Pomraning offered a new closure that is based on the time-independent asymptotic particles distribution due to the exact Boltzmann equation in an infinite, source-free homogeneous medium, which can be solved analytically~\cite{pomraning1964generalized}. As a matter of fact, this is a generalization of the asymptotic diffusion approximation~\cite{case1953introduction,CaseZweifel1967} (using $N=1$), for a general $N$. Pomraning originally derived this approximation in its neutron form however, we will derive it in its radiative transfer form. The last approximated equation takes the following form:
\begin{equation}\label{pomraningtruncation}
        \sigma_tI_N(x) + \left(\frac{N}{2N+1}\right)\pdv{I_{N-1}(x)}{x} + \left(\frac{N + 1}{2N + 1}\right) \pdv{(\alpha_N I_{N-1}(x))}{x} =0
\end{equation}
where $\alpha_N$ is the term that attempts to retain the symmetry of the transport equation~\cite{pomraning2005equations}, and is defined as:
\begin{equation}\label{alphapomraning0}
 \alpha_N = \frac{\int_{-1}^{1} I_a(x,\mu) P_{N+1}(\mu) d\mu}{\int_{-1}^{1} I_a(x,\mu) P_{N-1}(\mu) d\mu}
\end{equation}
$I_a$ is the spatial asymptotic form of the specific intensity. As one can notice, by setting $\alpha_N$ to zero, we obtain the classical $P_N$ approximation. In the general case, $I_a$ can be derived in multiple ways, each suitable for some angular distribution. Specifically, the asymptotic case can be derived from the exact solution of the time-independent, source-free Boltzmann equation for infinite homogeneous medium~\cite{case1953introduction,CaseZweifel1967}:
\begin{equation}\label{boltzmman_indep}
    \mu \pdv{ I(x,\mu)}{x} = \sigma_t \left(\frac{{\omega_{\mathrm{eff}}} }{2}cE(x) - I(x,\mu)\right)
\end{equation}
Where ${\omega_{\mathrm{eff}}}$ is the mean number of particles that are emitted per collision or source terms and is defined by (generalized to time-dependent case)~\cite{pomraning2005equations,cohen2018discontinuous}:
\begin{equation}
    {\omega_{\mathrm{eff}}}(x,t) = \frac{\sigma_s E(x,t) + \sigma_a B(x,t) + S(x,t)/c}{\sigma_t E(x,t)}.
\label{omega}
\end{equation}
${\omega_{\mathrm{eff}}}$ is the radiative transfer equivalence of the albedo $c$ that was used in~\cite{harel2020time} (not to be confused with the speed of light, in this paper), that is common mostly for neutronic modeling. Unlike in  our previous work~\cite{harel2020time}, in the case of strong sources ${\omega_{\mathrm{eff}}}$ may be much larger than unity (as an in multiplying media), which forces additional derivations for this work. Also notice that in the case of  LTE, ${\omega_{\mathrm{eff}}}\equiv 1$,  the asymptotic distribution tends to the diffusion limit.

By a separation of variables~\cite{case1953introduction,CaseZweifel1967}, one can solve the transport equation~\ref{boltzmman_indep} and obtain a closed transcendental equation for the eigenvalues of the positional asymptotic solution - $\varkappa_0$, as a function of ${\omega_{\mathrm{eff}}}$:
\begin{equation}
    \frac{2}{{\omega_{\mathrm{eff}}}}=\frac{1}{\varkappa_0}\ell n \left(\frac{1+\varkappa_0}{1-\varkappa_0}
    \right)
\label{kappa}
\end{equation}
and the asymptotic eigenfunctions - $I_a$ are (after forcing normalization):
\begin{equation}\label{asymptotc_eq}
I(x,\mu)\approx I_a(x,\mu)=\frac{{\omega_{\mathrm{eff}}}}{2}A_0\frac{e^{\varkappa_0\sigma_tx}}{1+\mu\varkappa_0}+
\frac{{\omega_{\mathrm{eff}}}}{2}B_0\frac{e^{-\varkappa_0\sigma_tx}}{1-\mu\varkappa_0}
\end{equation}

Plugging the asymptotic solution (Eq.~\ref{asymptotc_eq}) into the closure term (Eq.~\ref{alphapomraning0}), yields:
\begin{equation}\label{alphapomraning}
\alpha_N = \frac{\int_{-1}^{1} \left(\frac{1}{1+\mu\varkappa_0}\right) P_{N+1}(\mu) d\mu}{\int_{-1}^{1} \left(\frac{1}{1+\mu\varkappa_0}\right) P_{N-1}(\mu) d\mu}
\end{equation}
Plugging Eq.~\ref{alphapomraning} back into Eq.~\ref{pomraningtruncation} yields Pomraning's time-independent asymptotic $P_N$ approximation.

\subsection{The Time-dependent Asymptotic $P_N$}
Finally, we derive the time-dependent analogy to Pomranin'g closure, using the time-dependant source free homogeneous medium Boltzmann equation. This method is a the generalization for Heizler's asymptotic $P_1$ approximation~\cite{heizler2010asymptotic}, to a general $N$. The derivation is shown here briefly, focusing on the new features (for further details see~\cite{harel2020time}).

In order to derive the time-dependent asymptotic $P_N$ we start with Pomraning's closure, but with a slight modification in which we take the variable $\alpha_n$ outside the space derivative:
\begin{equation}\label{atd_nequation}
        \frac{1}{c}\pdv{I_N(x, t)}{t} + \sigma_t I_N(x,t) + \left(\frac{N}{2N+1}\right)\pdv{I_{N-1}(x, t)}{x} + \alpha_N\left(\frac{N + 1}{2N + 1}\right) \pdv{ I_{N-1}(x, t)}{x} = 0.
\end{equation}
This modification reproduces the asymptotic diffusion Fick's law, in the time-independent case for $N=1$~\cite{harel2020time}.
Additionally, $\alpha_n$ will now be derived by using the \textbf{time-dependent} source-free Boltzmann equation asymptotic distribution. Solving this equation is done by first applying Laplace transformation to Eq.~\ref{transport_Eq}:
\begin{equation}\label{s_boltzmann}
    \left(\frac{s}{c} + \sigma_t\right)\hat{I}_{s}(x,\mu) + \mu \pdv{ \hat{I}_{s}(x,\mu)}{x} = \frac{\sigma_t c {\omega_{\mathrm{eff}}} }{2}\hat{E}_{s}(x).
\end{equation}
Defining the following $s$-dependent coefficients:
\begin{subequations}\label{Hnotation}
    \begin{equation}
        \hat{\sigma}^{s}_t = \sigma_t + \frac{s}{c}
    \end{equation}
    \begin{equation}
        {\hat{\omega}_{\mathrm{eff}}}{}_{,s} = \frac{\sigma_s}{\hat{\sigma}^{s}_t} = \frac{{\omega_{\mathrm{eff}}}}{1 + \frac{s}{c\sigma_t}}
    \end{equation}
\end{subequations}
and plugging them back to Eq~\ref{s_boltzmann} yields:
\begin{equation}\label{boltzmann_I}
    \hat{\sigma}^{s}_t \hat{I}_{s}(x,\mu) + \mu \pdv{ \hat{I}_{s}(x,\mu)}{x} = \frac{\hat{\sigma}^{s}_t c {\hat{\omega}_{\mathrm{eff}}}{}_{,s} }{2}\hat{E}_{s}(x)
\end{equation}
Eq.~\ref{boltzmann_I} can be solved spatially in exactly the same manner as Eq.~\ref{boltzmman_indep}, but, in the Laplace domain. Therefore, the asymptotic solution is exactly the same except for the modified $s$-dependent coefficients:
\begin{equation}\label{TD_asymptotc_eq}
\hat{I}_{s}{}_{,a}(x,\mu)\approx\frac{{\hat{\omega}_{\mathrm{eff}}}{}_{,s}}{2}A_0\frac{e^{\hat{\varkappa}_{s,0}\hat{\sigma}^{s}_tx}}{1+\mu\hat{\varkappa}_{s,0}}+\frac{{\hat{\omega}_{\mathrm{eff}}}{}_{,s}}{2}B_0\frac{e^{-\hat{\varkappa}_{s,0}\hat{\sigma}^{s}_tx}}{1-\mu\hat{\varkappa}_{s,0}}
\end{equation}
With $s$-dependent closed transcendental equation for the asymptotic eigenvalues:
\begin{equation}
\frac{2}{{\hat{\omega}_{\mathrm{eff}}}{}_{,s}}=\frac{1}{\hat{\varkappa}_{s,0}}\ell n\left(\frac{1+\hat{\varkappa}_{s,0}}{1-\hat{\varkappa}_{s,0}}\right)
\label{hkappa}
\end{equation}
Plugging the time-dependent asymptotic solution - Eq.~\ref{TD_asymptotc_eq} into Eq.~\ref{alphapomraning0} yields:
\begin{equation}\label{alpha_s}
    \alpha_N\to\hat{\alpha}_{s}{}_{,N} = \frac{\int_{-1}^{1} \hat{I}_{s}{}_{,a}(x,\mu) P_{N+1}(\mu) d\mu}{\int_{-1}^{1} \hat{I}_{s}{}_{,a}(x,\mu) P_{N-1}(\mu) d\mu}= \frac{\int_{-1}^{1} \left(\frac{1}{1 + \mu \hat{\varkappa}_{s,0}}\right) P_{N+1}(\mu) d\mu}{\int_{-1}^{1} \left(\frac{1}{1 + \mu \hat{\varkappa}_{s,0}}\right) P_{N-1}(\mu) d\mu}
\end{equation}

Next, applying the Laplace transform to Eq. \ref{atd_nequation} using the $s$-dependent coefficients $\hat{\sigma}^{s}_t$ and $\hat{\alpha}_{s}{}_{,N}$ yields:
\begin{equation}\label{atd_nequation_laplace1}
    \hat{\sigma}^{s}_t\hat{I}_{s}{}_{,N}(x) + \left(\frac{N + (N+1)\hat{\alpha}_{s}{}_{,N}}{2N+1}\right)\pdv{\hat{I}_{s}{}_{,N-1}(x)}{x} = 0
\end{equation}
This Fick's law-form equation, defines the $s$-dependent ``diffusion" coefficient:
\begin{equation}\label{asymptotic_pn_ficklaw}
        \hat{I}_{s}{}_{,N}(x) = -\hat{D}_{s,N}(x) \pdv{\hat{I}_{s}{}_{,N-1}(x)}{x} \quad\longrightarrow\quad \hat{D}_{s,N}(x) = \frac{N+(N+1)\hat{\alpha}_{s}{}_{,N}}{(2N+1)\hat{\sigma}_t^s}
\end{equation}

As a generalization of Heizler's asymptotic $P_1$ approximation~\cite{heizler2010asymptotic}, an approximate closure (last $N$ equation), is derived by using two new coefficients ${\cal A}_N({\omega_{\mathrm{eff}}})$ and ${\cal B}_N({\omega_{\mathrm{eff}}})$:
\begin{equation}\label{atd_nequation_laplace2}
    \frac{{\cal A}_N({\omega_{\mathrm{eff}}})}{c} \pdv{I_N(x, t)}{t} + {\cal B}_N({\omega_{\mathrm{eff}}}) \sigma_t I_N(x,t) + \pdv{I_{N-1}(x,t)}{x} = 0.
\end{equation}
Further, applying the Laplace transform to the above approximated closure equation yields:
\begin{equation}\label{atd_nequation_laplace3}
    \hat{I}_{s}{}_{,N}(x) = - \frac{c}{s{\cal A}_N({\omega_{\mathrm{eff}}}) + c\sigma_t{\cal B}_N({\omega_{\mathrm{eff}}})} \cdot \pdv{\hat{I}_{s}{}_{,N-1}(x)}{x}
\end{equation}
and by substituting the approximate Eq.~\ref{atd_nequation_laplace3} into the exact closure Eq.~\ref{atd_nequation_laplace1} yields:
\begin{align}
& \frac{(2N+1)(c\sigma_t+s)}{N+(N+1)\hat{\alpha}_{s}{}_{,N}\left(\hat{\varkappa}_{s,0}({\hat{\omega}_{\mathrm{eff}}}{}_{,s})\right)}=\\
&\frac{(2N+1)(c\sigma_t+s)}{N+(N+1)\hat{\alpha}_{s}{}_{,N}\left(\hat{D}_0({\hat{\omega}_{\mathrm{eff}}}{}_{,s})\right)}\approx c\sigma_t{\cal B}_N({\omega_{\mathrm{eff}}}) + {\cal A}_N({\omega_{\mathrm{eff}}})s + {\cal O}(s^2)\nonumber.
\label{AnBn}
\end{align}
The modified asymptotic coefficients, ${\cal A}_N({\omega_{\mathrm{eff}}})$ and ${\cal B}_N({\omega_{\mathrm{eff}}})$ can be solved with the integral notation for $\hat{\alpha}_{s}{}_{,N}$ (Eq.~\ref{alpha_s}) under small-$s$ restrictions by using the Taylor expansion. $\hat{\alpha}_{s}{}_{,N}$ is an explicit function of $\hat{\varkappa}_{s,0}({\omega_{\mathrm{eff}}})$ or alternatively, $\hat{D}_0({\hat{\omega}_{\mathrm{eff}}}{}_{,s})$, for each ${\omega_{\mathrm{eff}}}$. For different functional expression for $\hat{\varkappa}_{s,0}({\omega_{\mathrm{eff}}})$ and $\hat{D}_0({\hat{\omega}_{\mathrm{eff}}}{}_{,s})$ please see~\cite{case1953introduction} and  the appendix of~\cite{harel2020time}.

The details of the derivation for $N=1$, 2 and 3, and also the derivation of $\hat{\alpha}_{s}{}_{,N}$ (which is practically finding the general solution for the integral in Eq.~\ref{alpha_s}) can be found in~\cite{harel2020time} (replacing ${\omega_{\mathrm{eff}}}$ in $c$) and we will not repeat it here. However we do mention, that for fully LTE problem, i.e. ${\omega_{\mathrm{eff}}}=1$, we reproduce the asymptotic $P_1$ results, ${\cal A}_1(1)=3/5$ and ${\cal B}_1(1)=3$~\cite{heizler2010asymptotic, heizler2012sp}. For asymptotic $P_2$ the method produces ${\cal A}_2(1)=4/7$ and ${\cal B}_2(1)=5/2$ and for asymptotic $P_3$ we yield ${\cal A}_3(1)=5/9$ and ${\cal B}_3(1) =7/3$. It is interesting to note that ${\cal A}_N(1)\to1/2$ and ${\cal B}_N(1)\to2$, when $N\to\infty$.

In Fig.~\ref{fig:AB} we present the behavior of the coefficients ${\cal A}_N({\omega_{\mathrm{eff}}})$ and ${\cal B}_N({\omega_{\mathrm{eff}}})$ for different values of ${\omega_{\mathrm{eff}}}$ for the range of $0\leqslant{\omega_{\mathrm{eff}}} \leqslant 5$, here for the first time for values that are greater than 1. This regime is needed for strong sources, inner or outer, far from LTE regimes that are important in radiative transfer. It can be used also for neutronics in multiplying regions (in the presence of fission cross-sections).
\begin{figure}[htbp!]
\includegraphics*[width=7.5cm]{A_bigc.eps}
\includegraphics*[width=7.5cm]{B_bigc.eps}
\caption{The media-dependent coefficients ${\cal A}_N({\omega_{\mathrm{eff}}})$ and ${\cal B}_N({\omega_{\mathrm{eff}}})$ as function of the medium's albedo ${\omega_{\mathrm{eff}}}$ for the first 3 orders of the time-dependent asymptotic $P_N$ approximation ($N=1,2,3$).}
\label{fig:AB}
\end{figure}

Fig.~\ref{fig:AB} suggests some important behaviors: First, ${\cal B}_N$ always equals to the inverse of Pomraning's time-independent asymptotic $P_N$ ``diffusion-like" coefficient, ${\cal B}_N=1/D_N({\omega_{\mathrm{eff}}})$. This forces the solution of the time-dependent asymptotic $P_N$ to tend to the solution of the time-independent asymptotic $P_N$ of Pomraning, when $t\to\infty$. Second, for fully-LTE medium (${\omega_{\mathrm{eff}}}=1$), ${\cal B}_N(1)=3$ equals to the classic $P_N$ approximation coefficient for general $N$ but ${\cal A}_N(1)$ is different (as in Heizler's asymptotic $P_1$ has shown~\cite{heizler2012sp}). Third, for pure absorbing case (in the limit of cold non-LTE regions ${\omega_{\mathrm{eff}}}=0$), ${\cal A}_N(0)={\cal B}_N(0)=1$ for any general $N$, yielding the exact particle velocity. We note that for ${\omega_{\mathrm{eff}}}>1$, the value of ${\cal A}_N$ is stabilized at ${\cal A}_N\approx 1/2$ for a wide range of ${\omega_{\mathrm{eff}}}$.

It should be noted that Pomraning proposed an approximate suggestion in~\cite{pomraning2005equations} (Eq. 3.147), which we call the $P_{1/{\cal B}_N}$ approximation. In this proposal, the {\em time-independent} closure is taken, therefore, determining ${\cal B}_N$ to be exactly as in the time-dependent asymptotic $P_N$ approximation:
\begin{equation}\label{pom_B}
\frac{1}{{\cal B}_N}\equiv D_N = \frac{N+(N+1)\alpha_N}{(2N+1)}
\end{equation}
when $\alpha_N$ is determined by Eq.~\ref{alphapomraning} and by setting ${\cal A}_N=1$.
For the case of {\em{purely steady-state case}} ${\omega_{\mathrm{eff}}}=1$, this $P_{1/{\cal B}_N}$ approximation has even simpler form, and it is the $N$'s expansion of the $P_{1/3}$ approximation, which was offered by Olson~\cite{olson2012alternate}, ${\cal A}_N=1$ and ${\cal B}_N=(2N+1)/N$ (which is the coefficient of the classic $P_N$ approximation when setting $\alpha_N=0$ in Eq.~\ref{pom_B}). The $P_{1/{\cal B}_N}$ approximation yields the exact particle velocity by definition due to the ${\cal A}_N=1$ choice. When approaching to the steady state, and as $N\to\infty$,  the asymptotic $P_N$ and the $P_{1/{\cal B}_N}$ approximations yields closer results. This means that the time-dependent asymptotic $P_N$ provides the rigorous mathematical justification (due to asymptotic analysis derivation) for the $P_{1/{\cal B}_N}$ approximation.

The results of Fig.~\ref{fig:AB} will serve as an input for any benchmark or simulation. In each zone (which in radiative transfer is any  numerical cell) ${\omega_{\mathrm{eff}}}(\vec{r},t)$ is calculated (via Eq.~\ref{omega}) for each time-step, due to the temperature dependency of the opacity and the difference between the material and radiation energies. Then we set ${\cal A}_N(\vec{r},t)$ and ${\cal B}_N(\vec{r},t)$ for each numerical cell in a given time-step by solving the linear set of $P_N$ equations, or simplified $P_N$ ($SP_N$) in multi-dimensions~\cite{spn}.

\section{Su-Olson Benchmark}\label{test_benchmark}

In order to test our time-dependent asymptotic $P_N$ approximation in its radiative form, we used the well-known Su-Olson benchmark ~\cite{su1997analytical}. The Su-Olson benchmark is a non-equilibrium gray one-dimensional radiative transfer source problem, with a constant opacity in an infinite medium.
\begin{figure}[htbp!]
\includegraphics*[width=8cm]{linear_c1_1.eps}
\includegraphics*[width=8cm]{linear_c1_2.eps}
\includegraphics*[width=8cm]{linear_c1_3.eps}
\caption{The radiation energy for the case of $c_s=0$ in linear scale. The different classic $P_N$ approximations are presented in dashed curves, while the asymptotic $P_N$ are presented by the solid curves. The exact benchmark results are shown in the black circles, while the $S_N$ results using $N=64$ are presented by the black curves. We have added the discontinuous asymptotic $P_1$ results from~\cite{cohen2018discontinuous} (the magenta curves).}
\label{fig:ca1_linear}
\end{figure}
\begin{figure}[htbp!]
\includegraphics*[width=8cm]{log_c1_1.eps}
\includegraphics*[width=8cm]{log_c1_2.eps}
\includegraphics*[width=8cm]{log_c1_3.eps}
\caption{The radiation energy for $c_s=0$ in log scale.}
\label{fig:ca1_log}
\end{figure}

It is common to present the physical properties as dimensionless as done by Cohen et al.~\cite{cohen2018discontinuous}, i.e. in dimensionless space units $z\equiv\sigma_{t} x$ and dimensionless time $\uptau\equiv\epsilon c \sigma_{t} t$. The material heat capacity is defined as: $C_v=\alpha T^3$ and $\epsilon=4a/\alpha$. It is also convenient to define the ratio of the scattering cross section to the total cross section $c_s=\sigma_s/ \sigma_{t}$, since we use dimensionless position variable.
The radiation source for this benchmark is isotropic and constant for $z \leqslant 0.5$ and $\uptau \leqslant 10$:
\begin{equation}
S(z,\uptau)=\begin{cases}
      1, & \text{if}\ \uptau\leq 10, \quad z\leq 0.5 \\
      0, & \text{otherwise}
    \end{cases}
    \label{source}
\end{equation}
We compared both the time-dependent asymptotic $P_N$ approximation to the classical $P_N$ approximation and included the benchmark analytic solutions and $S_N$ solution with $N=64$.
The numerical simulations for both the classical and the asymptotic $P_N$ approximation, were executed in a fully implicit scheme for time, with a finite difference scheme for space in {\em{Matlab}} (we have used the \textit{sparse} option for accelerating the band-matrix inversion). We used a constant spatial resolution of $\Delta z\approx 1.25\cdot 10^{-3}$ while $\Delta \uptau$ is defined dynamically such that the energy intensity will not change in each cell (between time steps) more than 0.05\%.

The first case that we used to test our approximation is the $c_s=0$, i.e. purely absorbing media in the absence of scattering. It is important to note that as opposed to neutronics, purely absorbing media does {\em{not}} mean ${\omega_{\mathrm{eff}}}=0$. In radiative transfer ${\omega_{\mathrm{eff}}}$ is determined by the ratio of the material energy to the radiation energy due to the black-body radiation source (and the external source as well). So in reaching LTE when $E\approx aT_m^4$ for sufficiently long times, ${\omega_{\mathrm{eff}}}\to1$ even in the purely absorbing case, due to the absorption-emission process.
\begin{figure}[htbp!]
\includegraphics*[width=8cm]{Mat_log_1.eps}
\includegraphics*[width=8cm]{Mat_log_2.eps}
\includegraphics*[width=8cm]{Mat_log_3.eps}
\caption{The material energy for $c_s=0$ in log scale.}
\label{fig:mat_ca1_log}
\end{figure}

In the following figures, the radiation energy as a function of the dimensionless space, is presented, where Fig.~\ref{fig:ca1_linear} is in a linear scale and Fig.~\ref{fig:ca1_log} is in a logarithmic scale. The linear scale focuses in the radiation in the bulk regions, while using the logarithmic scale we can test the tails of the distribution. The results are shown for $t=1,3.16$ and 10 dimensionless times. The different classic $P_N$ approximations are presented in dashed curves, while the asymptotic $P_N$ are presented by the solid curves. The exact benchmark results are shown in the black circles, while the $S_N$ results using $N=64$ are presented by the black curves.

As one can see, the benchmark exact solution and the $S_{64}$ fits perfectly, both in the bulk and in the tails. In the bulk regions, near the source (Fig~\ref{fig:ca1_linear}), the difference between classic and asymptotic $P_N$ is small, for any given $N$.
The major difference is near the tails of the distribution, as can be seen in the logarithmic scale (Fig~\ref{fig:ca1_log}). Each of the classical $P_N$ solutions yields a heat front that is too slow while the asymptotic time-dependent solution results in a heat front that are slightly faster. However, it is clear that the results produced by the asymptotic time-dependent approximation are better than the classical approximation, yielding almost the correct particle velocity.

Interesting feature is that the $P_2$ (both asymptotic and classic) have a discontinuous {\em{jump}} at the source surface, in $z=0.5$. This result, of the discontinuity of $P_2$ is well-known, due to the perpendicular direction of even $P_N$'s~\cite{RulkoLarsen1993}. It is quite similar to the jump of the discontinuous asymptotic $P_1$ results (the magenta curves)~\cite{cohen2018discontinuous}, however, the discontinuous asymptotic $P_1$ yields better results than both classic or asymptotic $P_2$, especially in long times. We note that such a discontinuity in the energy density is non-physical, since the exact solution of the RTE is of course, continuous. However, the physical justification of using discontinuities on the boundary between two different media rises, as mentioned above, from the physics of neutron transport. Discontinuous asymptotic diffusion theory ($P_1$ approximation) yields the correct critical radii in reactor-reflector systems, while the classic $P_N$ converges slowly compared to the exact solution~\cite{doyas}. This means that the non-physical existence of discontinuity is negligible compared to the benefits of correctly describing the particle's distribution  far from the boundary.

In Fig.~\ref{fig:mat_ca1_log} the complimentary results of the material energy is presented. The results of the material energy is similar qualitatively to the radiation energy results, of course differs quantitatively. As with the radiation energy, our new approximation yields better results than the classical one.

\begin{figure}[htbp!]
\includegraphics*[width=8cm]{log_c05_1.eps}
\includegraphics*[width=8cm]{log_c05_2.eps}
\includegraphics*[width=8cm]{log_c05_3.eps}
\caption{The radiation energy for $c_s=0.5$ in log scale.}
\label{fig:ca05_log}
\end{figure}
We also tested our new approximation for the scattering-included term of $c_s=0.5$. The following Fig~\ref{fig:ca05_log} presents the results. Similar conclusions regarding $c_s=0$ case, is valid for this one, however, the asymptotic $P_N$ results in this case are a little bit worse than in the pure absorbing case (mostly in early times), near the front of the heat waves.

\section{Discussion}
In this study, we presented the time-dependent asymptotic $P_N$ approximation for radiative transfer. The asymptotic $P_N$ approximation is based on the $P_N$ approximation and Pomraning's time-independent asymptotic closure method~\cite{pomraning1964generalized}, and was presented in previous work for general transport problems~\cite{harel2020time}. In this study we have focused on the modifications that had to be derived for radiative transfer problems, and the non-LTE radiative transfer benchmark problem, for which we examined this approximation.

The asymptotic $P_N$ approximation was derived in the same manner as all $P_N$ approximations are derived - by expanding the specific intensity into Legendre polynomials. Afterwards we introduced a new closure equation that is similar to Pomraning's closure.
As a generalization of the asymptotic $P_1$, the asymptotic $P_N$ approximation uses two media-dependent coefficients ${\cal A}_N({\omega_{\mathrm{eff}}})$ and ${\cal B}_N({\omega_{\mathrm{eff}}})$ {\em{in the closure equation}} that rests on the asymptotic distribution of the exact transport solution for infinite homogeneous medium, both in space (using Case et al., methodology~\cite{case1953introduction,CaseZweifel1967}) and time (using the Laplace domain in small-$s$). The coefficients ${\cal A}_N({\omega_{\mathrm{eff}}})$ and ${\cal B}_N({\omega_{\mathrm{eff}}})$ were derived for values of ${\omega_{\mathrm{eff}}} \leqslant 5$, thus enabling radiative transfer source problems in which ${\omega_{\mathrm{eff}}}$ may be in some cases, larger than unity (which may be also efficient for multiplying media in neutronics).

We have tested the new time-dependent asymptotic $P_N$ approximation for the well known Su-Olson benchmark in two cases - pure absorbing case ($c_s=0$) and scattering-included case $c_s=0.5$. The results show that our approximation share similar results to the classic $P_N$ in the bulk (near source) regions, and yields better results than the classical $P_N$ for any given $N$ near the tails of the distributions (the heat fronts regions). One important conclusion emerge from the results: the classic and the asymptotic $P_N$ approximations can be used as ``bounds" for the exact transport distribution, at least in the tails regions. The classic $P_N$ converges to the exact solution from the {\em{slower}} side, while the asymptotic $P_N$ converges from the {\em{faster}} side (though, converges faster in $N$ the classic $P_N$).

An interesting property, both the time-dependent $P_2$ (both classic and asymptotic) and the discontinuous asymptotic $P_1$ includes a jump at $z=0.5$. Although such discontinuity itself in the radiation energy density is of course non-physical, there are  benefits for describing correctly the  distribution of the particles far from the boundary.

In the future, we plan to extend our work to multi-dimensional problems, via the simplified $P_N$ ($SP_N$) approximation~\cite{spn}. Specifically, we are interested to explore the supersonic Marshak wave experiments, when a heat wave propagates into a low-Z foam, coated with a high-Z envelope (e.g. Au)~\cite{prr}. Since the coefficients ${\cal A}_N({\omega_{\mathrm{eff}}})$ and ${\cal B}_N({\omega_{\mathrm{eff}}})$ depend on the properties of numerical cells (i.e.~material), via ${\omega_{\mathrm{eff}}}(t,\vec{r})$ (which is function of time and space) and free of gradient (spatial derivatives) or time derivative, the extension to multi-dimensions is straightforward by using the existent time-dependent $SP_N$ code, which for low-order $N$, should yield a sufficiently good accuracy (same should work for other one-dimensional geometries). The only modification that is needed is setting the coefficients of the last (closure) equation to ${\cal A}_N({\omega_{\mathrm{eff}}})$ and ${\cal B}_N({\omega_{\mathrm{eff}}})$, that were derived and presented in this study (for example in Fig.~\ref{fig:AB}), instead of the classic $P_N$ coefficients.

\begin{acknowledgments}
We acknowledge the support of the PAZY Foundation under Grant \textnumero~61139927.
\end{acknowledgments}
\setlength{\baselineskip}{12pt}
\bibliographystyle{unsrtnat}