%~Mouliné par MaN_auto v.0.37.1 (6604bb9b) 2025-11-07 11:56:34
\documentclass[CRMATH,Unicode,biblatex,published]{cedram}

\TopicFR{\'Equations aux dérivées partielles}
\TopicEN{Partial differential equations}

\addbibresource{CRMATH_Bronsard_20241077.bib}

\usepackage{mathtools}
\usepackage{subcaption}

\newcommand{\dif}{\mathop{}\!{\operatorfont{d}}}

\DeclareMathOperator{\Rem}{Rem}
\DeclareMathOperator{\id}{id}
\DeclareMathOperator{\Id}{id}

%%%---------------------------------------------------------------------------

%%% DELIMITERS

\DeclarePairedDelimiter{\parens}{\lparen}{\rparen}
\DeclarePairedDelimiter{\braces}{\{}{\}}
\DeclarePairedDelimiter{\abs}{\lvert}{\rvert}
\DeclarePairedDelimiter{\norm}{\lVert}{\rVert}
\DeclarePairedDelimiter{\bracks}{[}{]}

%%%---------------------------------------------------------------------------

%%% ENVIRONEMENTS

%%% Env. {step}
\newcounter{step}
\newenvironment{step}{\refstepcounter{step}\begin{proof}[{Step~\thestep}]}{\let\qed\relax\end{proof}}

%%%---------------------------------------------------------------------------

%%% VERTICAL BAR

%%% ''Evaluated at'' macro
\newcommand\restr[3][]{{				
 \left.\kern-\nulldelimiterspace 
 #2 								
 \vphantom{\big\vert} 				
 \right\vert_{#3}^{#1}
}}

%%%---------------------------------------------------------------------------

\graphicspath{{./figures/}}

\newcommand*{\mk}{\mkern -1mu}
\newcommand*{\Mk}{\mkern -2mu}
\newcommand*{\mK}{\mkern 1mu}
\newcommand*{\MK}{\mkern 2mu}

\hypersetup{urlcolor=purple, linkcolor=blue, citecolor=red}

\newcommand*{\relabel}{\renewcommand{\labelenumi}{(\theenumi)}}
\newcommand*{\romanenumi}{\renewcommand*{\theenumi}{\roman{enumi}}\relabel}
\newcommand*{\Romanenumi}{\renewcommand*{\theenumi}{\Roman{enumi}}\relabel}
\newcommand*{\alphenumi}{\renewcommand*{\theenumi}{\alph{enumi}}\relabel}
\newcommand*{\Alphenumi}{\renewcommand*{\theenumi}{\Alph{enumi}}\relabel}
\let\oldtilde\tilde
\renewcommand*{\tilde}[1]{\mathchoice{\widetilde{#1}}{\widetilde{#1}}{\oldtilde{#1}}{\oldtilde{#1}}}

\title{Structure of Saturn ring defects for two small spherical colloidal particles}
\alttitle{Structure des défauts de type anneau de Saturne autour de deux petites particules colloïdales sphériques}

\author{\firstname{Lia} \lastname{Bronsard}\CDRorcid{0000-0002-0356-1714}}
\address{Department of Mathematics and Statistics, McMaster University, Hamilton, ON L8S 4L8 Canada}
\email{bronsard@mcmaster.ca}

\author{\firstname{Spencer} \lastname{Locke}\CDRorcid{0009-0002-7000-0447}}
\address{Department of Mathematics, University of Michigan, Ann Arbor, MI 48109, USA}
\email{lockespe@umich.edu}

\author{\firstname{Hayley} \lastname{Monson}}
\address[1]{Department of Mathematics and Statistics, McMaster University, Hamilton, ON L8S 4L8 Canada}
\email{monsonh@mcmaster.ca}

\author{\firstname{Dominik} \lastname{Stantejsky}\CDRorcid{0000-0001-7594-2922}}
\address{Universit\'e de Lorraine, Institut \'Elie Cartan de Lorraine, UMR 7502 CNRS, 54506 Vand\oe uvre-l\`es-Nancy Cedex, France}
\email{dominik.stantejsky@univ-lorraine.fr}

\subjclass{49K20, 35B38, 49S05}

\keywords{\kwd{Liquid crystals} \kwd{explicit solution} \kwd{line defect} \kwd{bispherical coordinates}}
\altkeywords{\kwd{Cristaux liquides} \kwd{solution explicite} \kwd{défaut de ligne} \kwd{coordonnées bisphériques}}

\thanks{L.B.\ was supported through NSERC Discovery Grant. H.M.\ and S.L.\ received funding from an NSERC-USRA grant. The main part of this work was carried out when S.L.\ and D.S.\ were affiliated with McMaster University}

\begin{abstract}
We study the small particle limit of the Landau--de Gennes model around two spherical colloids with strong homeotropic anchoring. We obtain an explicit representation of the minimizing $Q$-tensor. We then investigate the structure of defect lines and its dependence on the particle distance and orientation. In particular, for certain orientations and for small distances, we observe a line singularity disconnected from the singular line surrounding both particles, similar to the entangled hyperbolic defect configuration observed in experiments, while for larger distances, the two lines merge and eventually the inner singular line disappears.
\end{abstract}

\begin{altabstract}
Nous étudions la limite à petite échelle du modèle de Landau--de Gennes autour de deux colloïdes sphériques avec un ancrage homéotrope fort. Nous obtenons une représentation explicite du $Q$-tenseur minimisant. Nous étudions ensuite la structure des lignes de défaut et sa dépendance par rapport à la distance et \`a l'orientation entre les particules. En particulier, pour certaines orientations et pour de petites distances, nous observons une ligne singulière déconnectée de la ligne singulière entourant les deux particules, similaire à celle des défauts hyperboliques enchevêtrés observée dans les expériences, alors que pour de plus grandes distances, les deux lignes fusionnent et finalement la ligne singulière intérieure disparaît.
\end{altabstract}

\dateposted{2025-12-12}
\begin{document}

\maketitle

\section{Introduction}

We are interested in mathematically obtaining the structures of line defects around colloidal particles in a nematic liquid crystal. In the case of a single particle, this question has been studied for a variety of boundary conditions and limiting regimes by several authors, see e.g.~\cite{ABL,ABGL,ACS2021,ACS2024}. In the presence of more than one particle, one expects the particles to mutually influence each other, so we focus here on the model case of two spherical particles and their interaction. In the case of small particles and strong homeotropic anchoring (see the boundary conditions in~\eqref{eq:main-bvp}), we are able to obtain an explicit representation of the energy minimizing configuration, and this in turn allows us to study the singular structures analytically. In the case of a single spherical particle, this analysis was carried out in~\cite{ABL}. The energy that we are considering on the exterior of two balls $B_1$ and $B_2$ is given by
\begin{equation}\label{eq:energy}
E(Q) = \int_{\mathbb{R}^3\setminus (B_1\cup B_2)} \frac12 \abs{\nabla Q}^2 \dif x,
\end{equation}
which is obtained as the one constant approximation of the Landau--de Gennes energy functional in the limit of the particles being small compared to the nematic exchange length, see~\cite{ABL,Gartland}. Here, the map $Q$ is a $Q$-tensor that takes values into the space of $3\times 3$ symmetric matrices with zero traces. We note that in the case of unequal elastic constants, it is much harder to obtain analytical results even for a single spherical particle, see~\cite{BCMMSSV}.

The problem we study in this article frequently appears in the physical literature~\cite{ravnik09,musevicetal06,poulinstarklubenskyweitz97} where experiments and numerical simulations are carried out. They observe several different configurations for the singular structures that appear due to the colloidal presence in the liquid crystal. Here we present the mathematically rigorous results that we obtain in a limiting regime corresponding to particles being small compared to the nematic exchange length. The system of PDE that corresponds to critical points of this energy with homeotropic anchoring is as follows:
\begin{equation}\label{eq:main-bvp}
\left\{
\begin{aligned}
	-\Delta Q = \ & \mathbf{0} \quad && \mathbb{R}^{3}\backslash\overline{B_1(C_{1})\cup B_2(C_{2})},
\\	\restr{Q^{i}}{\partial B_{i}} = \ & \nu_{i}\otimes\nu_{i}-\frac{\Id}{3} && \partial B_i(C_{i}),i=1,2,
\\	Q^{i} \to \ & Q_{\infty} && \norm[\big]{ (x,y,z)} \to\infty,
\end{aligned}
\right.
%\begin{cases}
%	-\Delta Q=\mathbf{0} & \mathbb{R}^{3}\backslash\overline{B_1(C_{1})\cup B_2(C_{2})},
%\\	\restr{Q^{i}}{\partial B_{i}}=\nu_{i}\otimes\nu_{i}-\frac{\Id}{3} & \partial B_i(C_{i}),i=1,2,
%\\	Q^{i}\to Q_{\infty} & \norm[\big]{ (x,y,z)} \to\infty,
%\end{cases}
\end{equation}
where $Q_\infty= n_\infty \otimes n_\infty - \frac{1}{3}\id$ taking $n_\infty = \sin(\gamma)e_x + \cos(\gamma)e_z$ with $\gamma \in [0,\pi/2]$ determines the alignment at infinity, and where for $i=1,2$, we denote by $B_i$ the {balls} of centers $C_i$ on the $x$-axis and of radius~$r_i$, with $\nu_i$ their respective normals. We moreover define $D$ to be the distance between the two balls, so that $r_1 + r_2 + D$ gives the distance between the spheres' centers.

The solutions to this PDE that are obtained in Section~\ref{sec:expl-sol} can be visualized using ParaView~\cite{Paraview2015} as shown in Figures~\ref{fig:plot_pi2_D01-20}, \ref{fig:plot_pi2_3pi8_D05} and~\ref{fig:plot_0_3pi8_D05}, where the director fields and line defects are plotted. The director field is obtained from the $Q$-tensor by calculating the eigenvector to the largest eigenvalue, indicating the locally preferred axis of alignment of the liquid crystal molecules. When the two leading eigenvalues are equal, one loses the unique preferred direction and this is what yields a defect {in this model}. In our setting, we find singularities that are supported on lines of various geometries as depicted below, depending on the particle distance $D$, the radii $r_1$ and $r_2$, and the relative orientation of the particles with respect to $Q_\infty$.

\begin{figure}[h!]
\includegraphics[height=114pt]{2spheres_D01_pi2_n_angle.png}
~
\includegraphics[height=114pt]{2spheres_D20_pi2_n_angle.png}
\caption{Plot of the leading eigenvector of $Q$ (white lines) and the regions where the two leading eigenvalues are close (singularities, red lines). {Both balls have a radius~$1$ and the alignment at infinity is $n_\infty=n_z$ (i.e.\ $\gamma=0$)}. On the left, the distance $D$ between the balls is $D=0.1$ and there is one singular line surrounding both particles with one disconnected line in between them{, described in~\cite{Ravnik2007} as the \emph{entangled hyperbolic defect}}. On the right, $D=2$ and each particle is encircled by its own ring defect.}\label{fig:plot_pi2_D01-20}
\end{figure}

\begin{figure}[h!]
\includegraphics[height=114pt]{2spheres_D05_pi2_n_angle.png}
~
\includegraphics[height=114pt]{2spheres_D05_uneq_pi2_n_angle.png}
\caption{Plot of the leading eigenvector of $Q$ (white lines) and the regions where the two leading eigenvalues are close (singularities, red lines). The alignment at infinity has been chosen to be $n_\infty=n_z$ (i.e.\ $\gamma=0$). In both cases, the distance {$D$ between the balls is} $D=0.5$ and the line surrounding both particles meets the line in the center at two points. On the left, both particles have radius~$1$, while on the right, one particle has radius~$0.5$.}\label{fig:plot_pi2_3pi8_D05}
\end{figure}

\begin{figure}[h!]
\includegraphics[height=114pt]{2spheres_D05_0_n_angle.png}
~
\includegraphics[height=114pt]{2spheres_D05_3pi8_n_angle.png}
\caption{Plot of the leading eigenvector of $Q$ (white lines) and the regions where the two leading eigenvalues are close (singularities, red lines). Both balls have radius~$1$ and distance $D=0.5$. On the left, the particles are {positioned} along the axis {of alignment at infinity, i.e.\ $\gamma=\pi/2$}, while on the right, the two axes of alignment form an angle of $\frac{5}{8}\pi$. In both cases, each particle is encircled by its own ring defect, see also Figure~\ref{fig:plot_pi2_D01-20}.}\label{fig:plot_0_3pi8_D05}
\end{figure}

The main result of our analysis {is} for two particles of equal radius~$r_1=r_2=1$ aligned along the $x$-axis {with $\gamma=0$ and} is summarized in the numerical plot shown in Figure~\ref{fig:singular-struc}.

\begin{figure}
\includegraphics[width=0.8\linewidth]{Five-pts-diagram}
\caption{Plot showing the largest eigenvalues of $Q$ when restricted to the positive $y$-axis, as the distance $D$ varies. {The three pictures labeled~\textbf{A}, \textbf{B} and~\textbf{C} on the right visualize the three distinct configurations that can appear, depending on the distance $D$.}}\label{fig:singular-struc}
\end{figure}

As we can see, Figure~\ref{fig:singular-struc} shows numerically that for a given choice of colloidal distance $D$, there will be either zero, one, or two crossings of the leading eigenvalues of $Q$ along $y>0$, each of which corresponds to the globally observed singular structures shown in~\textbf{A}, \textbf{B}, and~\textbf{C} respectively. The structure~\textbf{A} is obtained for sufficiently large distances $D > D_c$, with $D_c$ indicated in Figure~\ref{fig:singular-struc} and estimated in Theorem~\ref{theo:singular-structure}, given that in this regime there are no leading eigenvalue crossings. This indicates the global defects become a pair of disjoint rings around the colloids for such values of~$D$. The structure~\textbf{B} is attained for distances slightly smaller than $D_c$, where there is a regime with a single crossing of the leading eigenvalues, thus demonstrating the presence of an inner ring intersecting the outer ring in the global picture of the defects. Finally, for small enough distances, the structure~\textbf{C} is achieved, given the presence of the regime where there are exactly two crossings of the leading eigenvalues, indicating that the inner ring has become disjoint from the outer ring in the global picture. Hence, this numerically demonstrates that the observed singular structure of the minimizer depends on the distance of the particles.

To prove the existence of these three regimes mathematically, we identify five points of this diagram that are stated and proved rigorously as points~\eqref{item:1}--\eqref{item:5} in Theorem~\ref{theo:singular-structure}, and which are labeled visually in Figure~\ref{fig:singular-struc}. To be specific, point~\eqref{item:1} confirms the conditions for the alignment at infinity, while point~\eqref{item:2} shows that for large distances, the Saturn rings around the two particles are separate, with no inner ring. Moreover, point~\eqref{item:3} implies that for small enough distances, there is an inner ring defect, and in fact that there is a critical distance denoted by $D_c$ such that, for any larger distances, no inner ring defects form.
Point~\eqref{item:4} shows that there is an open set of distances~$D$ below $D_c$ for which the inner ring intersects the outer ring. Finally, point~\eqref{item:5} demonstrates the existence of an additional critical distance, where for small enough distances, the inner ring will be disjoint from the outer ring.

The other elements of Figure~\ref{fig:singular-struc} are obtained numerically by approximating the series expansions for the eigenvalues in this regime (see equations~\eqref{eq:Q11-yaxis}, \eqref{eq:Q22-yaxis}, and~\eqref{eq:Q33-yaxis}) and analyzing the regimes where the largest ones cross. 
In Section~\ref{sec:bispherical} we introduce bispherical coordinates, which will allow us to represent the solutions of the system of Laplace equations on our domain. The solving process is carried out in Section~\ref{sec:expl-sol}, leading to an explicit series representation of each component. The series are then further analyzed in Section~\ref{sec:singularities} using Sturm's method for root finding in order to locate the defects as the crossing of principal eigenvalues.

\section{Bispherical coordinates and general solution to~(\ref{eq:main-bvp})}\label{sec:bispherical}

We introduce bispherical coordinates which are appropriate for the domain on which we study the energy. The bispherical coordinates are obtained via the rotation of a two-dimensional bipolar coordinate system about a focal line, creating a three dimensional orthogonal curvilinear coordinate system, see~\cite[p.~110]{FieldHandbook}.

The system corresponds to the coordinates $(\eta, \chi, \theta)$, made relative to the universal size parameter $0<a<\infty$ which is determined uniquely as a function of both radii and distance $D$. The coordinates define a unique point in $\mathbb{R}^3$, and are such that $-\infty < \eta < \infty$ identifies a sphere centered along the $x$-axis, while $\chi, \theta$ fix a point on the sphere, with $\chi \in [0,\pi]$ giving a point in the $xz$-plane, and $\theta \in [0,2\pi)$ giving the angle of rotation of the point about the $x$-axis. The coordinate transformation is defined by:
\[
x=\frac{a\sinh(\eta)}{\cosh(\eta)-\cos(\chi)},
\qquad y=\frac{a\sin(\chi)\sin(\theta)}{\cosh(\eta)-\cos(\chi)},
\qquad z=\frac{a\sin(\chi)\cos(\theta)}{\cosh(\eta)-\cos(\chi)},
\]
where the size parameter $a$ is such that there are two foci occurring at $x=\pm a$, acting as poles of the coordinate system, in the sense that $\eta \to \pm \infty$ corresponds to the spheres tending towards $x = \pm a$. In the context of this problem, we shall define the two spherical colloids to be balls $B_1$ and $B_2$ which have the radii $r_1$ and $r_2$, and that are centered at $C_1>0$ and $C_2 <0$ respectively on the $x$-axis.

It is shown in~\cite{CDH} that the size parameter $a$, and values $\eta_1, \eta_2$, are uniquely determined nonlinear functions of $r_1, r_2$, and $D$, such that $\eta = \eta_1$ and $\eta = \eta_2$ restricts to the surface of $B_1$ and $B_2$ respectively. Note that the region outside the two spheres is $\eta \in (\eta_2, \eta_1)$. In particular, as we shall utilize in Section~\ref{sec:singularities}, when fixing $r_1, r_2$ to be constant with the value $R>0$, the values of $a$, $\eta_1, \eta_2$ depend only on $D$. In fact, in this case, we can obtain $\eta_i = \eta_i(D)$ as a nonlinear function of $D$, which is given explicitly by~\cite{CDH} as:
\begin{equation}\label{eq:eta_1_exp}
\begin{cases}
	\eta_1 = \ln\parens[\Big]{\frac{2R+D + \sqrt{4RD + D^2}}{2R}},
\\	\eta_2 = -\eta_1.
\end{cases}
\end{equation}
Since the Laplacian is partially separable in the bispherical coordinates, the harmonic function $f \colon \mathbb{R}^3 \to \mathbb{R}$ solving $\Delta f = 0$ has the general solution given by
\begin{multline}\label{eq:bispherical-gen-sol}
f(\eta,\chi,\theta) =\parens[\big]{\cosh(\eta)-\cos(\chi)}^{\frac{1}{2}}\sum_{n=0}^{\infty}\sum_{m=0}^{n}\parens[\Big]{A_{m,n}\sinh\parens[\big]{\parens[\big]{n+\tfrac{1}{2}}\eta}+B_{m,n}\cosh\parens[\big]{\parens[\big]{n+\tfrac{1}{2}}\eta}}
\\	\cdot P_{n}^{m}\parens[\big]{\cos(\chi)}\parens[\big]{D_{m}\sin(m\theta)+E_{m}\cos(m\theta)},
\end{multline}
where $P_n^m$ are the associated Legendre polynomials; see~\cite[p.~111]{FieldHandbook}. The coefficients $A_{m,n},B_{m,n},D_m,E_n$ will be determined by the boundary conditions for our problem (see Section~\ref{sec:expl-sol}).

It is convenient to define the parameters $t_{\eta_1} \coloneqq e^{-\eta_1}$ and $t_{\eta_2} \coloneqq e^{\eta_2}$. In both cases, it can be verified that
\begin{equation}\label{eq:t-cosh-relation}
1-2t_{\eta_{i}}\cos(\chi)+t_{\eta_{i}}^{2}=2t_{\eta_{i}}\parens[\big]{\cosh(\eta_{i})-\cos(\chi)},
\end{equation}
for $i=1,2$. Note that since $\eta_i = \eta_i(D)$ by~\eqref{eq:eta_1_exp} when we have fixed both radii, we may also write $t_{\eta_i} = t_{\eta_i}(D)$ to be a function of $D$ when both radii are fixed. Moreover, taking the radii to be fixed, \eqref{eq:eta_1_exp} shows that $\eta_1 \to \infty$ and $\eta_2 \to -\infty$ as $D \to \infty$, as well $\frac{\dif \eta_1}{\dif D} > 0$ and $\frac{\dif \eta_2}{\dif D} < 0$ so that $\frac{\dif t_{\eta_i}}{\dif D} < 0$. We will also define the function $Z_n \colon \mathbb{R}^2 \to \mathbb{R}$ by
\begin{equation}\label{eq:Zn}
Z_n(\xi, \zeta) = \sinh\parens[\big]{\parens[\big]{n+\tfrac{1}{2}} \xi } - \cosh\parens[\big]{\parens[\big]{n+\tfrac{1}{2}}\xi}\tanh\parens[\big]{\parens[\big]{n+\tfrac{1}{2}}\zeta}.
\end{equation}
In order to calculate the explicit coefficients in~\eqref{eq:bispherical-gen-sol}, and thus obtain solutions to~\eqref{eq:main-bvp}, we will need the relations~\eqref{eq:gen-1}--\eqref{eq:gen-4}, \eqref{eq:gen-6} which follow by expanding the factors $\parens[\big]{\cosh(\eta) - \cos(\chi)}^{-1/2}$,  and its derivatives, in the required bases of associated Legendre polynomials. The associated Legendre polynomials possess the well-known (see~\cite[p.~773]{MathMethods}) generating function
\begin{equation}\label{eq:generating}
\sum_{n=0}^{\infty}q^{n}P_{n+m}^{m}\parens[\big]{\cos(\chi)}=\frac{(-1)^{m}(2m)! \, \sin^{m}(\chi)}{2^{m}m! \, \parens[\big]{1-2q\cos(\chi)+q^{2}}^{m+1/2}},
\end{equation}
for $\abs{q}<1$ that will be used to express the boundary conditions in terms of associated Legendre polynomials.

Observe that taking $m=0$, $m=1$, and $m=2$ in~\eqref{eq:generating}, with $q=t_{\eta_{i}}$ and recalling~\eqref{eq:t-cosh-relation} we have
\begin{align}
	\frac{1}{\parens[\big]{\cosh(\eta_{i})-\cos(\chi)}^{1/2}}&=\sum_{n=0}^{\infty}\sqrt{2}t_{\eta_{i}}^{n+1/2}P_{n}\parens[\big]{\cos(\chi)}, \label{eq:gen-1}
\\	\frac{\sin(\chi)}{\parens[\big]{\cosh(\eta_{i})-\cos(\chi)}^{3/2}} &=\sum_{n=1}^{\infty}2^{3/2}t_{\eta_{i}}^{n+1/2}P_{n}^{1}\parens[\big]{\cos(\chi)}, \label{eq:gen-2}
\\	\frac{\sin^2(\chi)}{\parens[\big]{\cosh(\eta_{i})-\cos(\chi)}^{5/2}} &=\sum_{n=2}^{\infty}\frac{2^{5/2}t_{\eta_{i}}^{n+1/2}}{3}P_{n}^{2}\parens[\big]{\cos(\chi)}, \label{eq:gen-3}
\end{align}
respectively. Moreover, taking $m=1$ in~\eqref{eq:generating}, differentiating with respect to $q$, multiplying by $2q/3$ and then adding this to the $m=1$ solution in~\eqref{eq:generating}, before evaluating everything at $q=t_{\eta_{i}}$ and employing~\eqref{eq:t-cosh-relation} reveals that, upon re-indexing of the sums
\begin{equation}\label{eq:gen-4}
\frac{\sin(\chi)}{\parens[\big]{\cosh(\eta_{i})-\cos(\chi)}^{5/2}}=\sum_{n=1}^{\infty}\frac{2^{5/2}(2n+1)t_{\eta_{i}}^{n+3/2}}{3(1-t_{\eta_{i}}^{2})}P_{n}^{1}\parens[\big]{\cos(\chi)}.
\end{equation}
To derive the final generating function we will need, we first apply the same steps used to derive~\eqref{eq:gen-4} except take $m=0$ and multiply by $2q$ instead of $2q/3$ to obtain the intermediary formula
\begin{equation}\label{eq:gen-5}
\frac{1}{\parens[\big]{\cosh(\eta_{i})-\cos(\chi)}^{3/2}}=\sum_{n=0}^{\infty}\frac{2^{3/2}(2n+1)t_{\eta_i}^{n+3/2}}{1-t_{\eta_{i}}^{2}}P_{n}\parens[\big]{\cos(\chi)}.
\end{equation}
Differentiating both sides of~\eqref{eq:gen-5}, inserting the formula $\frac{\dif}{\dif \chi}P_{n}\parens[\big]{\cos(\chi)}=\frac{n(n+1)}{(2n+1)\sin(\chi)}\parens[\big]{P_{n+1}\parens[\big]{\cos(\chi)}-P_{n-1}\parens[\big]{\cos(\chi)}}$ (see~\cite[p.~750]{MathMethods}) and re-indexing the summation shows that
\begin{equation}\label{eq:gen-6}
\frac{\sin^2(\chi)}{\parens[\big]{\cosh(\eta_{i})-\cos(\chi)}^{5/2}}=\frac{2^{5/2}}{3(1-t_{\eta_{i}}^{2})}\sum_{n=1}^{\infty}\parens[\big]{t_{\eta_{i}}^{n+5/2}(n+2)(n+1)-t_{\eta_{i}}^{n+1/2}n(n-1)}P_{n}\parens[\big]{\cos(\chi)}.
\end{equation}

\section{Derivation of explicit solutions}\label{sec:expl-sol}

In this section, we derive the explicit solutions for the minimizing $Q$-tensor by matching the coefficients $A_{m,n},B_{m,n},D_m,E_m$ in~\eqref{eq:bispherical-gen-sol} with the boundary conditions and the behavior at infinity.

Before calculating the solution to the PDE~\eqref{eq:main-bvp}, we first need two lemmas providing us with building blocks to construct the global solution using linearity of the system.

\begin{lemma}\label{lem:bvp-g_i}
The solution to the following boundary value problems for $i,j=1,2$ and $i\neq j$
\begin{equation}\label{eq:g_i-BVP}
\left\{
\begin{aligned}
	-\Delta g_{i} =	\ & 0	\quad	&& \mathbb{R}^{3}\backslash\overline{B_1(C_{1})\cup B_2(C_{2})},
\\	g_{i} =			\ & 1			&& \partial B_{i}(C_i),
\\	g_{i} =			\ & 0			&& \partial B_{j}(C_j),
\\	g_{i} \to		\ & 0			&& \norm[\big]{ (x,y,z)} \to\infty
\end{aligned}
\right.
%\begin{cases}
%	-\Delta g_{i} =	0	& \mathbb{R}^{3}\backslash\overline{B_1(C_{1})\cup B_2(C_{2})},
%\\	g_{i} = 1			& \partial B_{i}(C_i),
%\\	g_{i} = 0			& \partial B_{j}(C_j),
%\\	g_{i} \to 0			& \norm[\big]{ (x,y,z)} \to\infty
%\end{cases}
\end{equation}
is given explicitly by
\begin{equation}\label{eq:g_i-sol}
g_{i}(\eta,\theta,\chi)=\sqrt{2}\parens[\big]{\cosh(\eta)-\cos(\chi)}^{\frac{1}{2}}\sum_{n=0}^{\infty}t_{\eta_{i}}^{n+1/2}\frac{Z_{n}(\eta,\eta_{j})}{Z_{n}(\eta_{i},\eta_{j})}P_{n}\parens[\big]{\cos(\chi)}.
\end{equation}
\end{lemma}

\begin{proof}

The general solution for $g_{i}$ is given by~\eqref{eq:bispherical-gen-sol}. Since the boundary conditions have no $\theta$ dependence, we can choose $E_{m}=1$ for $m=0$, while we must have $D_{m}=E_{m}=0$ for all other~$m$. Hence, the general solution becomes
\[
g_{i}(\eta,\chi,\theta)=\sum_{n=0}^{\infty}\parens[\Big]{A_{n}\sinh\parens[\big]{\parens[\big]{n+\tfrac{1}{2}}\eta}+B_{n}\cosh\parens[\big]{\parens[\big]{n+\tfrac{1}{2}}\eta}}P_{n}\parens[\big]{\cos(\chi)}.
\]
The boundary condition $\restr{g_{i}}{\partial B_{i}}=g_{i}(\eta_{i},\chi,\theta)=1$ thus reads as
\[
\sum_{n=0}^{\infty}\parens[\Big]{A_{n}\sinh\parens[\big]{\parens[\big]{n+\tfrac{1}{2}}\eta_{i}}+B_{n}\cosh\parens[\big]{\parens[\big]{n+\tfrac{1}{2}}\eta_{i}}}P_{n}\parens[\big]{\cos(\chi)}=\frac{1}{\parens[\big]{\cosh(\eta_{i})-\cos(\chi)}^{1/2}}.
\]
Combining this with~\eqref{eq:gen-1} implies that, for all $n\geq0$
\begin{equation}\label{eq:g_i-sys-1}
A_{n}\sinh\parens[\big]{\parens[\big]{n+\tfrac{1}{2}}\eta_{i}}+B_{n}\cosh\parens[\big]{\parens[\big]{n+\tfrac{1}{2}}\eta_{i}}=\sqrt{2}t_{\eta_{i}}^{n+1/2}.
\end{equation}
Moreover, the boundary condition that $\restr{g_{i}}{\partial B_{j}}=g_{i}(\eta_{j})=0$ gives that
\begin{equation}\label{eq:zero-cond}
A_{n}\sinh\parens[\big]{\parens[\big]{n+\tfrac{1}{2}}\eta_{j}}+B_{n}\cosh\parens[\big]{\parens[\big]{n+\tfrac{1}{2}}\eta_{j}}=0.
\end{equation}
Solving the system of equations~\eqref{eq:g_i-sys-1} and~\eqref{eq:zero-cond} yields
\[
A_{n}=\frac{\sqrt{2}t_{\eta_{i}}^{n+1/2}}{Z_{n}(\eta_{i},\eta_{j})}
\qquad \text{and} \qquad
B_{n}=-A_n \tanh\parens[\big]{\parens[\big]{n+\tfrac{1}{2}}\eta_j}
\]
leading to~\eqref{eq:g_i-sol}.
\end{proof}

Next, we will solve a system related to~\eqref{eq:main-bvp} for $Q$. As $Q$ is a symmetric $3 \times 3$ matrix, we will solve component-wise for $Q_{\alpha \beta}^i$ where $i,j = 1,2$, $i\neq j$ and $\alpha, \beta = 1,2,3$ in the following lemma.

\begin{lemma}\label{lem:bvp-Q_r}
The solution to the following boundary value problem {for $i,j=1,2$ and $i\neq j$}
\begin{equation}\label{eq:Q_r-BVP}
\left\{
\begin{aligned}
	-\Delta Q^{i} = 	\ & \mathbf{0}							\quad	&& \mathbb{R}^{3}\backslash\overline{B_1(C_{1})\cup B_2(C_{2})},
\\	Q^{i} =			\ & \nu_{i}\otimes\nu_{i}-\frac{\Id}{3}			&& \partial B_{i}(C_i),
\\	Q^{i} =			\ & \mathbf{0}									&& \partial B_{j}(C_j),
\\	Q^{i} \to		\ & \mathbf{0} 									&& \norm[\big]{ (x,y,z)} \to\infty,
\end{aligned}
\right.
\end{equation}
are obtained component-wise. The off-diagonal components $Q_{12}^{i}$, $Q_{13}^{i}$, $Q_{23}^{i}$ of the solution are given explicitly by
\begin{equation}\label{eq:Q^i-off-diag}
Q_{\alpha\beta}^{i} = \parens[\big]{\cosh(\eta)-\cos(\chi)}^{1/2}\sum_{n=m_{\alpha\beta}}^{\infty}A_{n,i}^{\alpha\beta}\frac{Z_{n}(\eta,\eta_{j})}{Z_{n}(\eta_{i},\eta_{j})}P_{n}^{m_{\alpha\beta}}\parens[\big]{\cos(\chi)}f_{\alpha\beta}(\theta),
\end{equation}
where $f_{\alpha\beta}(\theta)=\sin(\theta)$ and $m_{\alpha\beta}=1$ for $Q_{12}^i$, $f_{\alpha\beta}(\theta)=\cos(\theta)$ and $m_{\alpha\beta}=1$ for $Q_{13}^i$, $f_{\alpha\beta}(\theta)=\sin(2\theta)$ and $m_{\alpha\beta}=2$ for $Q_{23}^i$, and where
\[
A_{n,i}^{12} = A_{n,i}^{13} =
\frac{a^{2}\sinh(\eta_{i})2^{5/2}(2n+1)t_{\eta_{i}}^{n+3/2}}{3r_{i}^{2}(1-t_{\eta_{i}}^{2})}-\frac{aC_{i}2^{3/2}t_{\eta_{i}}^{n+1/2}}{r_{i}^{2}}
\qquad \text{and} \qquad
A_{n,i}^{23} =
\frac{2^{3/2}a^{2}t_{\eta_{i}}^{n+1/2}}{3r_{i}^{2}}.
\]
Likewise, the diagonal components of the solution are given in the form
\begin{multline}\label{eq:Q^i-diag}
Q_{\alpha\alpha}^{i}=\parens[\big]{\cosh(\eta)-\cos(\chi)}^{1/2}
\\	\parens*{\sum_{n=0}^{\infty}A_{n,i}^{\alpha\alpha}\frac{Z_{n}(\eta,\eta_{j})}{Z_{n}(\eta_{i},\eta_{j})}P_{n}\parens[\big]{\cos(\chi)}+(-1)^{\alpha+1}\sum_{n=2}^{\infty}B_{n,i}^{\alpha\alpha}\frac{Z_{n}(\eta,\eta_{j})}{Z_{n}(\eta_{i},\eta_{j})}P_{n}^{2}\parens[\big]{\cos(\chi)}\cos(2\theta)},
\end{multline}
for $\alpha = 2,3$, where
\begin{align*}
	A_{n,i}^{33} & =A_{n,i}^{22}=\frac{2^{3/2}a^{2}}{3r_{i}^{2}(1-t_{\eta_{i}}^{2})}\parens[\big]{t_{\eta_{i}}^{n+5/2}(n+2)(n+1)-t_{\eta_{i}}^{n+1/2}n(n-1)}-\frac{\sqrt{2}}{3}t_{\eta_{i}}^{n+1/2},
\\	B_{n,i}^{11} & =B_{n,i}^{22}=A_{n,i}^{23}.
\end{align*}
The remaining components are found using the symmetry and tracelessness of the matrix, i.e.\ $Q_{21}^{i}=Q_{12}^{i}$, $Q_{31}^{i}=Q_{13}^{i}$, $Q_{32}^{i}=Q_{23}^{i}$, and $Q_{11}^i = -Q_{33}^i - Q_{22}^i$.
\end{lemma}

\begin{proof}
We generalize the method of generating functions used in the proof of Lemma~\ref{lem:bvp-g_i} to be able to handle the boundary condition
\[
\restr{Q^{i}}{\partial B_{i}}=\nu_{i}\otimes\nu_{i}-\frac{\Id}{3}=
\setlength\arraycolsep{5pt}
\begin{pmatrix}
	\nu_{x_{i}}^{2}-\tfrac{1}{3} 	& \nu_{x_{i}}\nu_{y_{i}} 		& \nu_{x_{i}}\nu_{z_{i}}
\\	\nu_{x_{i}}\nu_{y_{i}} 			& \nu_{y_{i}}^{2}-\tfrac{1}{3} 	& \nu_{y_{i}}\nu_{z_{i}}
\\	\nu_{x_{i}}\nu_{z_{i}} 			& \nu_{y_{i}}\nu_{z_{i}} 		& \nu_{z_{i}}^{2}-\tfrac{1}{3}
\end{pmatrix},
\]
where
\[
\nu=(\nu_{x_{i}},\nu_{y_{i}},\nu_{z_{i}})=\parens[\Bigg]{\frac{1}{r_{i}}\parens[\Bigg]{\frac{a\sinh(\eta_{i})}{\cosh(\eta_{i})-\cos(\chi)}-C_{i}},\frac{a\sin(\chi)\sin(\theta)}{r_{i}\parens[\big]{\cosh(\eta_{i})-\cos(\chi)}},\frac{a\sin(\chi)\cos(\theta)}{r_{i}\parens[\big]{\cosh(\eta_{i})-\cos(\chi)}}}
\]
is the unit normal to the spheres in the bispherical coordinates.

\paragraph{Solution for \boldmath{$Q_{12}^i$}}
Using the general solution~\eqref{eq:bispherical-gen-sol}, and that the only dependence on $\theta$ in the boundary conditions for $Q_{12}^i$ is $\sin(\theta)$, we see that $D_{m}=1$ for $m=1$, while $D_{m}=E_{m}=0$ for all other $m$. Thus,
\[
Q_{12}^{i}(\eta,\chi,\theta)=\parens[\big]{\cosh(\eta)-\cos(\chi)}^{1/2}\sum_{n=1}^{\infty}\parens[\Big]{A_{n}\cosh\parens[\big]{\parens[\big]{n+\tfrac{1}{2}}\eta}+B_{n}\sinh\parens[\big]{\parens[\big]{n+\tfrac{1}{2}}\eta}}P_{n}^{1}\parens[\big]{\cos(\chi)}\sin(\theta).
\]
Using~\eqref{eq:gen-2} and~\eqref{eq:gen-4} with the condition that $\restr{Q_{12}^{i}}{\partial B_{i}}=\nu_{x_{i}}\nu_{y_{i}}$ implies that
\begin{equation}\label{Q_12-sys-1}
A_{n}\cosh\parens[\big]{\parens[\big]{n+\tfrac{1}{2}}\eta_{i}}+B_{n}\sinh\parens[\big]{\parens[\big]{n+\tfrac{1}{2}}\eta_{i}}=\frac{2^{5/2}a^{2}\sinh(\eta_{i})(2n+1)t_{\eta_{i}}^{n+3/2}}{3r_{i}^{2}(1-t_{\eta_{i}}^{2})}-\frac{aC_{i}2^{3/2}t_{\eta_{i}}^{n+1/2}}{r_{i}^{2}}.
\end{equation}
Since $\restr{Q_{12}^{i}}{\partial B_{j}}=0$ yields~\eqref{eq:zero-cond}, combining this with~\eqref{Q_12-sys-1}, we obtain the coefficients
\[
A_{n}=\frac{2^{5/2}a^{2}\sinh(\eta_{i})(2n+1)t_{\eta_{i}}^{n+3/2}}{3r_{i}^{2}(1-t_{\eta_{i}}^{2})Z_n(\eta_i,\eta_j)}-\frac{aC_{i}2^{3/2}t_{\eta_{i}}^{n+1/2}}{r_{i}^{2}Z_n(\eta_i,\eta_j)}
\qquad \text{and} \qquad
B_{n}=-A_n \tanh\parens[\big]{\parens[\big]{n+\tfrac{1}{2}}\eta_j},
\]
thus yielding~\eqref{eq:Q^i-off-diag} in the case $\alpha = 1, \beta = 2$.

\paragraph{Solution for \boldmath{$Q_{13}^{i}$}}
Since $Q_{13}^{i}$ solves the same equation as $Q_{12}^{i}$, with the only difference in the boundary condition $\restr{Q_{13}^{i}}{\partial B_{i}}=\nu_{x_{i}}\nu_{z_{i}}$ being that the factor of $\sin(\theta)$ in $\restr{Q_{12}^{i}}{\partial B_{i}}=\nu_{x_{i}}\nu_{y_{i}}$ becomes a factor of $\cos(\theta)$, the same procedure shows that the solution to $Q_{13}^{i}$ is the same as $Q_{12}^{i}$, except with this change of factors.

\paragraph{Solution for \boldmath{$Q_{23}^i$}}
Using the general solution~\eqref{eq:bispherical-gen-sol}, and that the only dependence on $\theta$ in the boundary conditions for $Q_{23}^i$ is $\sin(\theta)\cos(\theta) =\frac12 \sin(2\theta)$, we must necessarily have that $D_{m}=1$ for $m=2$, while $D_{m}=E_{m}=0$ for all other $m$. Thus,
\[
Q_{23}^{i}(\eta,\chi,\theta)=\parens[\big]{\cosh(\eta)-\cos(\chi)}^{1/2}\sum_{n=2}^{\infty}\parens[\Big]{A_{n}\sinh\parens[\big]{\parens[\big]{n+\tfrac{1}{2}}\eta}+B_{n}\cosh\parens[\big]{\parens[\big]{n+\tfrac{1}{2}}\eta}}P_{n}^{2}\parens[\big]{\cos(\chi)}\sin(2\theta).
\]
Applying~\eqref{eq:gen-3} to the condition that $\restr{Q_{23}^{i}}{\partial B_{i}}=\nu_{y_{i}}\nu_{z_{i}}$, and using $\sin(\theta)\cos(\theta) =\frac12 \sin(2\theta)$, reveals that
\begin{equation}\label{eq:Q_23-sys-1}
A_{n}\sinh\parens[\big]{\parens[\big]{n+\tfrac{1}{2}}\eta_{i}}+B_{n}\cosh\parens[\big]{\parens[\big]{n+\tfrac{1}{2}}\eta_{i}}=\frac{2^{3/2}a^{2}t_{\eta_{i}}^{n+1/2}}{3r_{i}^{2}}.
\end{equation}
The equation $\restr{Q_{23}^{i}}{\partial B_{j}}=0$ yields the same condition as~\eqref{eq:zero-cond}, and thus combining this with~\eqref{eq:Q_23-sys-1} we obtain the coefficients
\begin{equation}
A_{n}=\frac{2^{3/2}a^{2}t_{\eta_{i}}^{n+1/2}}{3r_{i}^{2}Z_{n}(\eta_{i},\eta_{j})}
\qquad \text{and} \qquad
B_{n}=-A_n \tanh\parens[\big]{\parens[\big]{n+\tfrac{1}{2}}\eta_j},
\end{equation}
yielding~\eqref{eq:Q^i-off-diag} in the case $\alpha = 2$, $\beta = 3$.

\paragraph{Solution for \boldmath{$Q_{33}^{i}$}}
Using the general solution~\eqref{eq:bispherical-gen-sol}, and noting that since the boundary conditions
\begin{equation}\label{eq:Q_33_BC}
\restr{Q_{33}^{i}}{\partial B_{i}}=\nu_{z_{i}}^{2}-\frac{1}{3}=\parens[\Bigg]{\frac{a^{2}\sin^2(\chi)}{2r_{i}^{2}\parens[\big]{\cosh(\eta_{i})-\cos(\chi)}^{2}}-\frac{1}{3}}+\frac{a^{2}\sin^2(\chi)\cos(2\theta)}{2r_{i}^{2}\parens[\big]{\cosh(\eta_{i})-\cos(\chi)}^{2}}
\end{equation}
involve a term with a $\cos(2\theta)$ factor and a term independent of $\theta$, we see that $E_{m}=1$ for $m=0,2$, while $D_{m}=E_{m}=0$ for all other $m$. Thus,
\begin{multline*}
Q_{33}^{i}(\eta,\chi,\theta) =\parens[\big]{\cosh(\eta)-\cos(\chi)}^{1/2}\parens[\Bigg]{\sum_{n=0}^{\infty}\parens[\Big]{A_{n}^{0}\cosh\parens[\big]{\parens[\big]{n+\tfrac{1}{2}}\eta}+B_{n}^{0}\sinh\parens[\big]{\parens[\big]{n+\tfrac{1}{2}}\eta}}P_{n}\parens[\big]{\cos(\chi)}
\\	+\sum_{n=2}^{\infty}\parens[\Big]{A_{n}^{2}\cosh\parens[\big]{\parens[\big]{n+\tfrac{1}{2}}\eta}+B_{n}^{2}\sinh\parens[\big]{\parens[\big]{n+\tfrac{1}{2}}\eta}}P_{n}^{2}\parens[\big]{\cos(\chi)}\cos(2\theta)}.
\end{multline*}
Now, using~\eqref{eq:gen-1} and~\eqref{eq:gen-6} for the first term in~\eqref{eq:Q_33_BC}, and~\eqref{eq:gen-3} for the second term in~\eqref{eq:Q_33_BC} shows
\begin{multline*}
A_{n}^{0}\cosh\parens[\big]{\parens[\big]{n+\tfrac{1}{2}}\eta_{i}} + B_{n}^{0}\sinh\parens[\big]{\parens[\big]{n+\tfrac{1}{2}}\eta_{i}}
\\	=	\frac{2^{3/2}a^{2}}{3r_{i}^{2}(1-t_{\eta_{i}}^{2})}\parens[\big]{t_{\eta_{i}}^{n+5/2}(n+2)(n+1)-t_{\eta_{i}}^{n+1/2}n(n-1)}-\frac{\sqrt{2}}{3}t_{\eta_{i}}^{n+1/2},
\end{multline*}
and
\[
A_{n}^{2}\cosh\parens[\big]{\parens[\big]{n+\tfrac{1}{2}}\eta_{i}}+B_{n}^{2}\sinh\parens[\big]{\parens[\big]{n+\tfrac{1}{2}}\eta_{i}}=\frac{a^{2}2^{3/2}t_{\eta_{i}}^{n+1/2}}{3r_{i}^{2}}
\]
for $n\geq2$. The boundary condition $\restr{Q_{33}^{i}}{\partial B_{i}}=0$ implies that~\eqref{eq:zero-cond} holds for the coefficients $A_{n}^{0},B_{n}^{0}$ as well as $A_{n}^{2},B_{n}^{2}$. Solving the respective systems yields the coefficients satisfying~\eqref{eq:Q^i-diag} with $\alpha = 3$.

\paragraph{Solution for \boldmath{$Q_{22}^{i}$}}
Finally, since
\[
\restr{Q_{22}^{i}}{\partial B_{i}}=\nu_{y_{i}}^{2}-\frac{1}{3}=\parens[\Bigg]{\frac{a^{2}\sin^2(\chi)}{2r_{i}^{2}\parens[\big]{\cosh(\eta_{i})-\cos(\chi)}^{2}}-\frac{1}{3}}-\frac{a^{2}\sin^2(\chi)\cos(2\theta)}{2r_{i}^{2}\parens[\big]{\cosh(\eta_{i})-\cos(\chi)}^{2}},
\]
it follows that $Q_{22}^{i}$ solves the same equation as $Q_{33}^{i}$, with the same boundary condition as $\restr{Q_{33}^{i}}{\partial B_i}$ except with the factor of $\cos(2\theta)$ changed to $-\cos(2\theta)$. Thus, the same procedure shows that the solution to $Q_{22}^{i}$ is the same as $Q_{33}^{i}$, except with this change of factors.
\end{proof}

\begin{rema}\label{rema:new-coeffs}
It follows by construction of the bispherical coordinate system, and is verifiable using the exact formulas for $a, \eta_1, \eta_2$ in~\cite{CDH}, that $a^{2}/r_{i}^{2}=\sinh^2(\eta_{i})$. Hence, one can write all the coefficients in terms of the parameters $\eta_i$, and obtain for example
\[
A_{n,i}^{23}=\frac{(t_{\eta_{i}}^{2}+t_{\eta_{i}}^{-2}+2)t_{\eta_{i}}^{n+1/2}}{3\sqrt{2}},A_{n,i}^{33}=\frac{1-t_{\eta_{i}}^{2}}{3\sqrt{2}t_{\eta_{i}}^{2}}\parens[\big]{t_{\eta_{i}}^{n+5/2}(n+2)(n+1)-t_{\eta_{i}}^{n+1/2}n(n-1)}-\frac{\sqrt{2}}{3}t_{\eta_{i}}^{n+1/2}.
\]
The other coefficients can also be obtained in terms of the same parameters.
\end{rema}

Combining Lemmas~\ref{lem:bvp-g_i} and~\ref{lem:bvp-Q_r} proves the following theorem.

\begin{theo}\label{theo:total-Q}
The solution $Q$ of the Dirichlet boundary value problem
\begin{equation}\label{eq:total-PDE}
\left\{
\begin{aligned}
	-\Delta Q =					\ & \mathbf{0}	\quad 	&& \mathbb{R}^{3}\backslash\overline{B_1(C_{1})\cup B_2(C_{2})},
\\	\restr{Q^{i}}{\partial B_{i}} =	\ & \nu_{i}\otimes\nu_{i}-\frac{\Id}{3}	&& \partial B_i(C_{i}), \ i=1,2,
\\	Q \to						\ & Q_{\infty} 		&& \norm[\big]{(x,y,z)} \to\infty,
\end{aligned}
\right.
%\begin{cases}
%	-\Delta Q =	\mathbf{0}												& \mathbb{R}^{3}\backslash\overline{B_1(C_{1})\cup B_2(C_{2})},
%\\	\restr{Q^{i}}{\partial B_{i}} = \nu_{i}\otimes\nu_{i}-\frac{\Id}{3}	& \partial B_i(C_{i}), \ i=1,2,
%\\	Q \to Q_{\infty} 													& \norm[\big]{(x,y,z)} \to\infty,
%\end{cases}
\end{equation}
is given by $Q=Q^{1}+Q^{2}+(1-g_{1}-g_{2})\, Q_{\infty}$, {where $g_{i}$ solves~\eqref{eq:g_i-BVP}, $Q^{i}$ solves~\eqref{eq:Q_r-BVP} and $Q_\infty$ is given by
\[
Q_{\infty}=n_{\infty}\otimes n_{\infty}-\frac{1}{3}\id,
\]
for $n_{\infty}=\sin(\gamma)e_{x}+\cos(\gamma)e_{z}$\ and $\gamma \in [0,\pi/2]$.}
\end{theo}

\begin{nota*}
For a convergent series of functions $S(x) = \sum_{n=m}^\infty s_n(x)$, for $N\geq m$ we shall use the notation $S_N(x)$ to denote the first $N$ terms in the sum, and $\Rem(S,N)$ to denote the remainder such that
\begin{equation}\label{eq:series-decomp}
S(x) = S_N(x) + \Rem\parens[\big]{S(x),N}.
\end{equation}
\end{nota*}

\begin{rema}[Asymptotic bounds]\label{rema:asymptotics}
One can obtain general bounds on the remainders of each of the components in equations~\eqref{eq:Q^i-off-diag} and~\eqref{eq:Q^i-diag}. Indeed, since $\eta \in (\eta_2, \eta_1)$, then $\abs[\big]{Z_{n}(\eta,\eta_{j})/$ $Z_{n}(\eta_{i},\eta_{j})}\leq1$. Moreover one can bound the Legendre polynomials with $\abs[\big]{P_{n}(\cos x)}\leq K(\chi)^2$, for $K(\chi)=\parens[\big]{1+(\pi^{4}/16)\sin^{2}(\chi)}^{-1/16}$, see~\cite{ElbertLaforgia}, and $\abs[\big]{P_{n}^{1}(\cos x)}\leq2\pi^{-1/2}(n+1)nK(\chi)$, $\abs[\big]{P_{n}^{2}(\cos x)}\leq2\pi^{-1/2}(n+1)n(n-1)K(\chi)$, which were shown in~\cite{Lohofer}. Therefore, it follows by applying the triangle inequality, along with these estimates, that we may bound the remainders of the off-diagonal components \eqref{eq:Q^i-off-diag} as
\[
\abs[\big]{\Rem(Q_{\alpha\beta}^{i},N)}\leq2\pi^{-1/2}\parens[\big]{\cosh(\eta)-\cos(\chi)}^{1/2}\abs[\big]{f_{\alpha\beta}(\theta)}K(\chi)R_{\alpha\beta}^{N}(t_{\eta_{i}}),
\]
while we may bound the remainders of the diagonal components \eqref{eq:Q^i-diag} for $\alpha=2,3$ with
\begin{equation}\label{eq:bound-diag}
\abs[\big]{\Rem(Q_{\alpha\alpha}^{i},N)}\leq\parens[\big]{\cosh(\eta)-\cos(\chi)}^{1/2}K(\chi)\parens[\Big]{K(\chi) R_{\alpha\alpha}^{N}(t_{\eta_{i}})+2\pi^{-1/2}\abs[\big]{\cos(2\theta)}R_{23}^{N}(t_{\eta_{i}})},
\end{equation}
for some infinite series $R_{\alpha \beta}^N$, and where $f_{\alpha\beta}(\theta)$ is as defined in Lemma~\ref{lem:bvp-Q_r}. These infinite series can be evaluated explicitly as they can all be reduced to combinations of derivatives of the standard geometric series using that $\abs{t_{\eta_i}}<1$. In particular, these bounds are finite for all colloidal distances $D>0$. Finally, recalling $Q_{11}^{i}=-Q_{22}^{i}-Q_{33}^{i}$, one can use~\eqref{eq:bound-diag} to obtain a bound on $Q_{11}^{i}$. In the proof of Theorem~\ref{theo:singular-structure} in Section~\ref{sec:singularities}, we derive explicit forms for $R_{\alpha \beta}^N$, which we use to obtain precise error estimates.
\end{rema}

\begin{rema}[Equal radii] \label{rema:sol-equal-radii}
The explicit solutions of the minimizing $Q$-tensor simplify dramatically in the case where the colloids have equal radii. Given the symmetry of the system, this corresponds to the case where $\eta_{1}=-\eta_{2}$, so that by~\eqref{eq:Zn}, $Z_{n}(\eta_{i},\eta_{j})=2\sinh\parens[\big]{\parens[\big]{n+\tfrac{1}{2}}\eta_{i}}$. Hence, it follows by Theorem~\ref{theo:total-Q} that in the case of equal radii for $\alpha,\beta\in\braces{1,2,3} $ with $\alpha\neq\beta$,
\begin{equation}\label{eq:eq-radii-off-diag}
Q_{\alpha\beta}=\parens[\big]{\cosh(\eta)-\cos(\chi)}^{1/2}\sum_{n=m_{\alpha\beta}}^{\infty}A_{n}^{\alpha\beta}\frac{\cosh\parens[\big]{\parens[\big]{n+\tfrac{1}{2}}\eta}}{\cosh\parens[\big]{\parens[\big]{n+\tfrac{1}{2}}\eta_{1}}}P_{n}^{m_{\alpha\beta}}\parens[\big]{\cos(\chi)}f_{\alpha\beta}(\theta)
,
\end{equation}
where $f_{\alpha\beta}(\theta)=\sin(\theta)$ and $m_{\alpha\beta}=1$ for $Q_{12}$, $f_{\alpha\beta}(\theta)=\cos(\theta)$ and $m_{\alpha\beta}=1$ for $Q_{13}$, $f_{\alpha\beta}(\theta)=\sin(2\theta)$ and $m_{\alpha\beta}=2$ for $Q_{23}$.
Here, we note that in the case of equal radii, $A_{n}^{\alpha\beta} \coloneqq A^{\alpha\beta}_{n,1} = A^{\alpha\beta}_{n,2}$ from Lemma~\ref{lem:bvp-Q_r}. Likewise, when $\alpha=\beta \neq 1$ the diagonal components become
\begin{multline}\label{eq:eq-radii-diag}
Q_{\alpha\alpha} =\parens[\big]{\cosh(\eta)-\cos(\chi)}^{1/2}\left(\sum_{n=0}^{\infty}A_{n}^{\alpha\alpha}\frac{\cosh\parens[\big]{\parens[\big]{n+\tfrac{1}{2}}\eta}}{\cosh\parens[\big]{\parens[\big]{n+\tfrac{1}{2}}\eta_{1}}}P_{n}\parens[\big]{\cos(\chi)}\right.
\\	\left. + (-1)^{\alpha+1}\sum_{n=2}^{\infty}B_{n}^{\alpha\alpha}\frac{\cosh\parens[\big]{\parens[\big]{n+\tfrac{1}{2}}\eta}}{\cosh\parens[\big]{\parens[\big]{n+\tfrac{1}{2}}\eta_{1}}}P_{n}^{2}\parens[\big]{\cos(\chi)}\cos(2\theta)\right)-\frac{1}{3}\parens{1-g_{1}-g_{2}}, 
\end{multline}
where $B^{\alpha\alpha}_n \coloneqq B^{\alpha\alpha}_{n,1} = B^{\alpha\alpha}_{n,2}$ are as defined in Lemma~\ref{lem:bvp-Q_r}, and $Q_{11} = -Q_{33} - Q_{22}$.
\end{rema}

\section{Singular structures for equal radii}\label{sec:singularities}

In this section, we prove Theorem~\ref{theo:singular-structure} in order to understand the singular structure of the solution for two equal radii and $n_\infty=e_z$, i.e.\ $\gamma = 0$. This theorem shows existence of  points~\eqref{item:1}--\eqref{item:5} from Figure~\ref{fig:singular-struc}, thereby justifying the configurations~\textbf{A}, \textbf{B} and~\textbf{C}. As Figure~\ref{fig:singular-struc} demonstrates, it is enough to restrict $Q$ to the $y$-axis to distinguish when the three observed configurations~\textbf{A}, \textbf{B} and~\textbf{C} occur. Hence, we shall restrict our analysis to this case.

In bispherical coordinates, the $y$-axis corresponds to $\eta = 0$, $\theta = \pi/2$. For simplicity, we place the colloids symmetrically about the origin along the $x$-axis, so that $C_1 = -C_2 > 0$. With the equal radii, we note that, for any fixed colloidal distance $D$, we have $t_{\eta_1} = t_{\eta_2}$ given that $t_{\eta_1} = e^{-\eta_1} = e^{\eta_2} = t_{\eta_2}$, and so we may simply drop the subscript to simplify the notation. Note that since we are working with the two radii $r_1,r_2$ fixed, the function $t$ introduced in Section~\ref{sec:bispherical} is a function solely of $D$. In fact, one may determine the explicit non-linear relation between $t$ and~$D$ in the case of $r_1 = r_2 = 1$ using the formula~\eqref{eq:eta_1_exp}, which reads
\begin{equation}\label{eq:t(D)}
t(D) = \frac{2}{(D+2) + \sqrt{D^2+4D}}.
\end{equation}
We also write $t$ instead of $t(D)$ to improve readability whenever possible.

As outlined in Remark~\ref{rema:sol-equal-radii}, the solution {components} are given as follows
\begin{align}
	& Q_{11}(\chi,t) =\parens[\big]{1-\cos(\chi)}^{1/2}\sum_{n=0}^{\infty}\frac{\parens{-2\tilde{A}_{n}+\sqrt{2}t^{n+1/2}}}{\cosh\parens[\big]{\parens[\big]{n+\tfrac{1}{2}}\eta_{1}}}P_{n}\parens[\big]{\cos(\chi)}-\frac{1}{3}, \label{eq:Q11-yaxis}
\\	&
\begin{multlined}[b][.9\displaywidth]\label{eq:Q22-yaxis}
Q_{22}(\chi,t) = \parens[\big]{1-\cos(\chi)}^{1/2}
\\	\parens*{\sum_{n=0}^{\infty}\frac{\tilde{A}_{n}}{\cosh\parens[\big]{\parens[\big]{n+\tfrac{1}{2}}\eta_{1}}}P_{n}\parens[\big]{\cos(\chi)}+\sum_{n=2}^{\infty}\frac{B_{n}^{33}}{\cosh\parens[\big]{\parens[\big]{n+\tfrac{1}{2}}\eta_{1}}}P_{n}^{2}\parens[\big]{\cos(\chi)}}-\frac{1}{3},
\end{multlined}
\\	&
\begin{multlined}[b][.9\displaywidth]\label{eq:Q33-yaxis}
Q_{33}(\chi,t) = \parens[\big]{1-\cos(\chi)}^{1/2}
\\	\parens*{\sum_{n=0}^{\infty}\frac{\parens{\tilde{A}_{n}-\sqrt{2}t^{n+1/2}}}{\cosh\parens[\big]{\parens[\big]{n+\tfrac{1}{2}}\eta_{1}}}P_{n}\parens[\big]{\cos(\chi)}-\sum_{n=2}^{\infty}\frac{B_{n}^{33}}{\cosh\parens[\big]{\parens[\big]{n+\tfrac{1}{2}}\eta_{1}}}P_{n}^{2}\parens[\big]{\cos(\chi)}}+\frac{2}{3},
\end{multlined}
\end{align}
where $\tilde{A}_{n}=A_{n}^{33}+\frac{\sqrt{2}}{3}t^{n+1/2}$, while the other components vanish.

Since in this case the solution $Q$ is diagonal, the eigenvectors are given by the canonical basis and the eigenvalues are just the diagonal entries. As the singularities occur precisely at the locations where the two leading eigenvalues are equal, it is enough to distinguish the regions in which each of the basis vectors are the dominating eigenvector.

\begin{theo}\label{theo:singular-structure}
For $r_{1}=r_{2}=1$, $\eta=0$ and $\theta=\pi/2$, $Q$ is diagonal and there exists $D_{c}>D_{2}>D_{1}>D_{m}>0$ such that:
\begin{enumerate}
\item \label{item:1}
for $y\to\infty$, the dominating eigenvector is $e_{z}$;
\item \label{item:2}
for $D\to\infty$, the dominating eigenvector on the $y$-axis is $e_{z}$;
\item \label{item:3}
for $D_{m}<D<D_{c}$, $e_{x}$ is dominating at $x=(0,0,0)$, while for $D>D_{c}$ $e_{z}$ is dominating at $x = (0,0,0)$;
\item \label{item:4}
for $D=D_{2}$, $e_{y}$ is never the dominating eigenvector on the $y$-axis;
\item \label{item:5}
for $D=D_{1}$, there exists a point $y_{1}>0$ such that $e_{y}$ is the dominating eigenvector at $x=(0,y_{1},0)$.
\end{enumerate}
{Furthermore,} $D_{c}\in(2.01,2.03)$, and one can take $D_{2}=1.85$, $D_{1}=0.35$ and $D_{m}=\frac{4}{\sqrt{3}}-2$.
\end{theo}

\begin{proof}
We will prove each point individually.
\begin{proof}[Proof of~(\ref{item:1})]
Taking $y\to\infty$ means $\chi\to 0$ in the bispherical coordinates. Thus, with $\chi
\to 0$ in~\eqref{eq:Q11-yaxis}, \eqref{eq:Q22-yaxis}, and~\eqref{eq:Q33-yaxis}, using the fact that $P_n(1)=1$, $P_n^2(1)=0$, we see that $Q_{11}, Q_{22} \to -1/3$, $Q_{33} \to 2/3$ as $\chi \to 0$, which completes the proof of~\eqref{item:1}.
\let\qed\relax
\end{proof}

\begin{proof}[Proof of~(\ref{item:2})]
Since in the bispherical system $D\to\infty$ corresponds to $\eta_{1}\to\infty$ and $t\to0$, it follows using~\eqref{eq:Q11-yaxis}, \eqref{eq:Q22-yaxis}, and~\eqref{eq:Q33-yaxis} as well as Remark~\ref{rema:new-coeffs} that $\tilde{A}_{n}\to0$ for $n\geq 0$ and $B_{n}^{33}\to0$ when $n\geq2$, and thus we see that $Q_{11}(\chi,t)\to-\frac{1}{3}$, $Q_{22}(\chi,t)\to-\frac{1}{3}$ and $Q_{33}(\chi,t)\to\frac{2}{3}$ as $D\to\infty$, which completes the proof of~\eqref{item:2}.
\let\qed\relax
\end{proof}

\begin{proof}[Proof of~(\ref{item:3})]
Taking $y=0$ corresponds to $\chi=\pi$ in bispherical coordinates, and hence $P_{n}\parens[\big]{\cos(\chi)}=(-1)^{n}$ and $P_{n}^{2}\parens[\big]{\cos(\chi)}=0$. It follows by~\eqref{eq:Q22-yaxis} and~\eqref{eq:Q33-yaxis} that $Q_{33}>Q_{22}$, so that $Q_{33}$ is dominant if and only if $Q_{33}>Q_{11}$. Thus, we will study the roots of $Q_{33}-Q_{11}$, as a function of $D$, which can be computed as
\begin{equation}\label{eq:S(D)}
S(t)\coloneqq(Q_{33}-Q_{11})(t)=1+\sum_{n=0}^{\infty}(-1)^{n} h(n,t),
\end{equation}
where {the function $h(n,t)$ is defined as}
\begin{equation}\label{eq:terms-in-sum-S}
h(n,t)\coloneqq\frac{2t^{n+1/2}(1-t^{2})\parens[\big]{t^{n+5/2}(n+2)(n+1)-t^{n+1/2}n(n-1)}-8t^{2n+3}}{t^{2}(1+t^{2n+1})},
\end{equation}
and where we used Remark~\ref{rema:new-coeffs} and the fact that $1/\cosh({N\eta_i}) = 2t^N / (t^{2N} + 1)$. Recall that $t$ is a function of $D$.

We will prove~\eqref{item:3} in several steps.

\begin{step}\label{step_1}
\emph{Showing $S(D)<0$ for $D\in[D_{m},1.96]$.}

Recalling the {decomposition}~\eqref{eq:series-decomp}, we proceed by bounding $S_3$ and $\Rem(S,3)$. First, for $S_3$, an explicit computation with~\eqref{eq:S(D)} and~\eqref{eq:terms-in-sum-S} reveals that
\[
\frac{\dif S_{3}}{\dif t}=\frac{L_{1}(t)}{(1-t+t^{2}+t^{5}-t^{6}+t^{7})^{2}},
\]
where $L_{1}(t)$ is a degree~15 polynomial in $t$. An application of Sturm's Theorem (see~\cite[p.~52]{Sturm}) shows that $L_{1}(t)$ has exactly two roots $t\in(0,1)$, with the smallest root being $t<0.6330$, with $L_{1}(t)>0$ for $t<0.6330$. Using~\eqref{eq:t(D)}, this means that $D>0.22$, since we are taking $D>D_{m}$ and we recall $D_m \coloneqq \frac{4}{\sqrt{3}}-2>0.22$. Thus, $S_{3}(D)$ is increasing for all $D>D_{m}$.

Next, we will bound $\Rem(S,3)$. First observe that the equation $h(n,t)=0$ implies
\[
(n+2)(n+1)t^{4}-2(n^{2}+n-1)t^{2}+n(n-1) = 0,
\]
whose smallest root $t\in(0,1)$ is $\tau_{*}=\sqrt{\frac{n-1}{n+1}}$ for $n\geq2$. One can directly verify that for $t=\frac{1}{\sqrt{n+1}}<\tau_{*}$, it holds $h(n,t)<0$. Therefore, it follows that for any $t\in(0,\tau_{*})$ and $n\geq2$ with $n$ even, $h(n,t)(-1)^{n}<0$ while $h(n,t)(-1)^{n}>0$ for $n$ odd. Using the bispherical formula~\eqref{eq:t(D)} it can be shown that $D>D_{m}$ implies $t(D)\in(0,\tau_{*})$.

Therefore, it holds
\begin{equation}\label{eq:upper-bound-Rem-S-4}
\begin{split}
\Rem(S,3)
	& \leq-\sum_{j=0}^{\infty}h\parens[\big]{3+2j,t(D)}
\\	& \leq-\sum_{j=0}^{\infty}2(1-t^2)\parens[\big]{t^{4j+7}(2j+5)(2j+4)-t^{4j+5}(3+2j)(2j+2)}-8t^{4j+7}
\\	& = \frac{L_{2}(t)}{(1-t^{2})^{2}(1+t^{2})^{3}} \eqqcolon U(t),
\end{split}
\end{equation}
where $L_{2}(t)$ is a degree~15 polynomial in $t$. Differentiating this upper bound yields
\[
U'(t)=\frac{\tilde{L}_{2}(t)}{(1-t^{2})^{2}(1+t^{2})^{4}}
\]
where $\tilde{L}_{2}(t)$ is a degree~16 polynomial. Applying Sturm's Theorem to this polynomial shows that it has no roots in the interval $(0,1)$, and thus since one can compute that $\tilde{L}_{2}(0.5)<0$, it follows that $U'(t)>0$ for all $t\in(0,1)$. Thus, recalling that $t=t(D)$ via~\eqref{eq:t(D)} and $\dif t/\dif D < 0$, the chain rule shows that $\frac{\dif}{\dif D}\bracks[\big]{U\parens[\big]{t(D)}}<0$, and thus $U\parens[\big]{t(D)}$ is monotonically decreasing in $D$.

We can now use the fact that $S_3$ is increasing and $U$ is decreasing in distance $D$ to show that $S\parens[\big]{t(D)} < 0$ for $D \in [D_m,1.96]$. Because monotonicity on the whole interval $[D_m,1.96]$ only indicates that $S\parens[\big]{t(D)} < S_3\parens[\big]{t(1.96)} + U\parens[\big]{t(D_m)}$, which is computed to be a positive quantity, then we will instead divide the interval into two sub-intervals and estimate $S\parens[\big]{t(D)}$ on both of them, finding that the monotonicity provides a negative upper bound.

First assuming that $D\in[D_{m},1.46]$, the monotonicity results show
\[
S\parens[\big]{t(D)}<S_{3}\parens[\big]{t(1.46)}+U\parens[\big]{t(D_m)}+1<0\,.
\]
Repeating this process, but assuming now $D\in[1.46,1.96]$, shows that $S\parens[\big]{t(D)}<0$ as well.
\end{step}

\begin{step}
\emph{Strict monotonicity of $S$ for $D\geq1.96$.}

We begin by noting that by chain rule
\[
\frac{\dif}{\dif D}S\parens[\big]{t(D)}=\frac{\dif t}{\dif D}\cdot\parens[\Bigg]{\frac{\dif S_{3}}{\dif t}+\frac{\dif}{\dif t}\Rem\parens[\Big]{S\parens[\big]{t(D)},3}},
\]
and we will bound each of the terms. First, since it can be shown that $\dif S_3/\dif t$ is not monotonic in the desired region, we will use Sturm's method directly to bound it. In particular, we will show $\dif S_3/\dif t < -3$ and thus we will look at the roots of the difference $\dif S_3/\dif t + 3$ to show it is strictly negative. By explicitly differentiating and rearranging, one can show that
\[
\frac{\dif S_3}{\dif t} + 3
\ = \
\frac{Q(t)}{(1+t)^2(1-t+t^2)^2(1-t+t^2-t^3+t^4)^2},
\]
where $Q(t)$ is a degree~15 polynomial in $t$. An application of Sturm's Theorem shows that $Q(t)<0$ for all $t\in (0,0.3)$, which implies that $Q\parens[\big]{t(D)}<0$ for all $D \geq 1.96$ by~\eqref{eq:t(D)}. Hence, $\frac{\dif S_3}{\dif t} \parens[\big]{t(D)}< -3$ for all $D \geq 1.96$. Moreover, by differentiating $\Rem(S,3)$ across the sum, and hence {differentiating}~\eqref{eq:terms-in-sum-S}, one may obtain the bound
\[
\begin{split}
\frac{\dif}{\dif t}\Rem(S,3)
	& < \sum_{n=3}^{\infty}2t^{2n-2}(12+2n^{3}+24n)
\\	& = \frac{4(-38t^{10}+133t^{8}-158t^{6}+69t^{4})}{(1-t^{2})^{4}}\eqqcolon U_{1}(t).
\end{split}
\]
With the same procedure that was applied to~\eqref{eq:upper-bound-Rem-S-4}, one can use Sturm's Theorem to demonstrate that $U_1\parens[\big]{t(D)}$ is monotonically decreasing in $D$, so $U_1\parens[\big]{t(D)} \leq U_1\parens[\big]{t(1.96)} < 2$ by exact computations. It follows that $\frac{\dif S_{3}}{\dif t}+\frac{\dif}{\dif t}\Rem(S,3)<2-3<0$, and hence since $\frac{\dif t}{\dif D}<0$ then $S$ is strictly monotonically increasing in $D$. Thus, there can exist at most one root.
\end{step}

\begin{step}
\emph{Bounding $D_c$.}

Using these facts, one can approach the value of $D_{c}$ by estimating $S\parens[\big]{t(D)}$ in nested intervals that are decreasing in length. Taking $D\in(1.96,2.01)$, we have
\[
S\parens[\big]{t(D)}<S_{3}\parens[\big]{t(2.01)}+U\parens[\big]{t(2.01)}<0
\]
by direct computation. Hence, combining this with Step~\ref{step_1} shows that $S\parens[\big]{t(D)} < 0$ for $D_m < D < 2.01$.

Moreover, taking $S\parens[\big]{t(D)}=S_{4}\parens[\big]{t(D)}+\Rem\parens[\big]{S\parens[\big]{t(D)},4}$, one can show that $S\parens[\big]{t(D)}>0$ for $D>2.03$. With the same method to obtain the upper bound for $\Rem(S,3)$ in~\eqref{eq:upper-bound-Rem-S-4}, one can show that
\[
\Rem(S,4) \geq\sum_{j=1}^{\infty}h(2+2j,t)\geq\frac{L_{3}(t)}{(1-t^{2})^{2}(1+t^{2})^{3}}\eqqcolon L(t)
\]
where $L_{3}(t)$ is a degree~17 polynomial that can be computed. Therefore, for $D>2.03$ monotonicity gives
\[
S\parens[\big]{t(D)}>S_{4}\parens[\big]{t(2.03)}+L\parens[\big]{t(2.03)}>0,
\]
where we used an explicit calculation of the functions to determine the final bound. Hence, there must exist a unique root $D_{c}\in(2.01,2.03)$, finishing the proof of~\eqref{item:3}.
\end{step}
\let\qed\relax
\end{proof}

\begin{proof}[Proof of~(\ref{item:4})]
Fixing $D=D_{2}=1.85$, we will show that $Q_{22}(\chi)-Q_{33}(\chi)<0$ for all $\chi\in[0,\pi]$. Let us observe that
\begin{equation}\label{eq:Q_22-Q_33}
\begin{split}
Q_{22}(\chi,t) - Q_{33}(\chi,t)
	& =\sum_{n=2}^{\infty}\frac{2E(\chi)B_{n}^{11}}{\cosh\parens[\big]{\parens[\big]{n+\tfrac{1}{2}}\eta_{1}}}P_{n}^{2}(\chi)+\sum_{n=0}^{\infty}\frac{\sqrt{2}E(\chi)t^{n+1/2}}{\cosh\parens[\big]{\parens[\big]{n+\tfrac{1}{2}}\eta_{1}}}P_{n}(\chi)-1
\\	& \eqqcolon f_{1}(\chi)+f_{2}(\chi)-1, 
\end{split}
\end{equation}
where $E(\chi)\coloneqq (1-\cos \chi)^{1/2}$. To complete the proof, we will obtain specific bounds on $f_1(\chi)$ and $f_2(\chi)$ to show that~\eqref{eq:Q_22-Q_33} is negative.

We begin by bounding $f_1(\chi)$. Using the bound $P_{n}^{2}\parens[\big]{\cos(\chi)}\leq2\pi^{-1/2}(n+2)(n+1)n(n-1)$ from~\cite{Lohofer}, one observes that
\begin{equation}\label{eq:Remainder-f}
\begin{split}
\Rem\parens[\big]{f_{1}(\chi),2}
	& \leq \frac{\sinh^{2}(\eta_{1})2^{9/2}E(\chi)}{3\sqrt{\pi}}\sum_{n=4}^{\infty}t^{2n+1}(n+2)(n+1)n(n-1)
\\	& = \frac{\sinh^{2}(\eta_{1})2^{9/2}E(\chi)}{3\sqrt{\pi}}\frac{L_{4}(t)}{(1-t^{2})^{5}}
\end{split}
\end{equation}
where $L_{4}(t)$ is a degree~17 polynomial in $t$. Using the fact that $E(\chi)\leq\sqrt{2}$ in~\eqref{eq:Remainder-f}, and evaluating the resulting expression exactly at $D=D_2$, we have $\Rem\parens[\big]{f_{1}(\chi),2}<0.08$ and hence
\[
f_{1}(\chi) \leq 2E(\chi)B_{2}^{11}P_{2}^{2}\parens[\big]{\cos(\chi)}+2E(\chi)B_{3}^{11}P_{3}^{2}\parens[\big]{\cos(\chi)}+0.08.
\]
Expanding the associated Legendre polynomials in this expression allows us to write this bound in the form $f_{1}(\chi) \leq E(\chi)h(\chi)+0.08$, where $h(\chi)$ is a degree~3 polynomial in $\cos(\chi)$ satisfying $h(\chi)<0.06$ for all $\chi\in[0,\pi]$, which can be verified by solving $h'(\chi) = 0$, a quadratic in $\cos(\chi)$, to obtain the maximum point. This leads to the bound $f_{1}(\chi)<0.16$. Next, we will bound $f_2(\chi)$. Performing a similar procedure as in~\eqref{eq:Remainder-f} for $f_{2}(\chi)$, we note that
\[
\Rem\parens[\big]{f_{2}(\chi),2} \leq \frac{4t^{5}}{1-t^{2}}<0.008,
\]
where we evaluated at $D=D_2$. Hence we now write $f_{2}(\chi)\leq E(\chi)k(\chi)+0.008$ where $k(\chi)$ is a linear function in $\cos(\chi)$, and the function $E(\chi)k(\chi)$ is maximized at $\chi=\pi$ with $E(\chi)k(\chi)<0.8$. Hence, it follows that $f_{1}(\chi)+f_{2}(\chi)<1$, which completes the proof of~\eqref{item:4}.
\let\qed\relax
\end{proof}

\begin{proof}[Proof of~(\ref{item:5})]
Observe that Figure~\ref{fig:singular-struc} shows that at $D = D_1$, there is a whole interval of $y$-values, and thus corresponding $\chi$-values, where $e_y$ is dominating. However, for simplicity, and to show existence, we will choose $\chi=\pi/3$. We must show that $(Q_{22} - Q_{33})(\pi/3, t) > 0$ and $(Q_{22} - Q_{11})(\pi/3, t) > 0$. Using~\eqref{eq:Q_22-Q_33}, we observe that $Q_{22}(\pi/3, t) - Q_{33}(\pi/3, t) = f_1(\pi/3) + f_2(\pi/3) - 1$. Using an identical bounding technique as in~\eqref{eq:Remainder-f}, except with the optimized bound $\abs[\big]{P_{n}^{2}(1/2)}\leq\sqrt{\frac{8}{\pi}}\parens[\big]{\frac{4}{3}}^{5/4}(n+2)(n+1)$ from~\cite{Lohofer}, one can see that
\begin{equation}\label{eq:rem-f1}
\abs[\big]{\Rem(f_{1}(\pi/3),5)}\leq\sqrt{\frac{8}{\pi}}\parens[\Bigg]{\frac{4}{3}}^{5/4}\frac{8\sinh^2(\eta_1)}{3}\frac{L_4(t)}{(1-t^{2})^{3}}<0.0166,
\end{equation}
where $L_4(t)$ is a degree~21 polynomial in $t$, and the last bound is obtained by evaluating the function at $D = D_1$, recalling that for simplicity we noted $\eta_1 = \eta_1(D)$ and $t = t(D)$. Hence, $f_1(\pi/3) \geq {(f_1)}_{N=5}(\pi/3) - 0.0166 > 0.20$ using an exact evaluation on the finite sum to obtain the bound. Note that with this same method used for bounding the series in~\eqref{eq:rem-f1}, but instead employing Bernstein's inequality $\abs[\big]{P_{n}(1/2)}\leq(2/\pi)^{1/2}(4/3)^{1/4}$ from~\cite{Lohofer}, one can show that $\abs[\big]{\Rem(f_{2}(\pi/3),5)}<0.005$ and thus $f_2(\pi/3) \geq {(f_2)}_{N=5}(\pi/3) - 0.005 > 0.829$. Therefore, we see that $f_{1}(\pi/3)+f_{2}(\pi/3)>1.029$ so that $Q_{22}(\pi/3,t)>Q_{33}(\pi/3,t)$.

Similar bounding techniques will show that $Q_{22}(\pi/3,t)>Q_{11}(\pi/3,t)$. In particular, we can express
\begin{equation}\label{eq:Q22-Q11_pi/3}
\begin{split}
Q_{22}(\pi/3,t)-Q_{11}(\pi/3,t)
	& =\sum_{n=0}^{\infty}\frac{3\tilde{A}_{n}}{\sqrt{2}\cosh\parens[\big]{\parens[\big]{n+\tfrac{1}{2}}\eta_{1}}}P_{n}\parens[\big]{\cos(\pi/3)}
\\	& \quad + \sum_{n=2}^{\infty}\frac{B_{n}^{11}}{\sqrt{2}\cosh\parens[\big]{\parens[\big]{n+\tfrac{1}{2}}\eta_{1}}}P_{n}^{2}\parens[\big]{\cos(\pi/3)}
\\	& \quad - \sum_{n=0}^{\infty}\frac{t^{n+1/2}}{\cosh\parens[\big]{\parens[\big]{n+\tfrac{1}{2}}\eta_{1}}}P_{n}\parens[\big]{\cos(\pi/3)}
\\	& \eqqcolon f_{3}(\pi/3)+\frac{1}{2}f_{1}(\pi/3)-f_{2}(\pi/3),
\end{split}
\end{equation}
where $f_{1},f_{2}$ were defined in~\eqref{eq:Q_22-Q_33} and $f_{3}$ refers to the series on the first line. We will demonstrate how one can bound the terms in this sum to show the quantity is positive.

Since we found $\abs[\big]{\Rem\parens[\big]{f_2(\pi/3),5}} < 0.005$, then we see that
\[
-f_{2}(\pi/3)\geq-{(f_2)}_{N=5}(\pi/3) -\abs[\big]{\Rem\parens[\big]{f_2(\pi/3),5}}>-0.835-0.005=-0.84.
\]
Likewise, we also have that $f_{1}(\pi/3)/2>0.1$, and hence all that remains is to bound $f_{3}(\pi/3)$. Again employing Bernstein's inequality, we obtain an analogous bound to~\eqref{eq:rem-f1} showing $\abs[\big]{\Rem(f_{3}(\pi/3),8)}<0.005$ so that
\[
f_{3}(\pi/3)\geq {(f_3)}_{N=8}(\pi/3)-0.005>0.746-0.005=0.741,
\]
where the bounds are computed with the exact evaluation of the functions. Thus, inserting these bounds into~\eqref{eq:Q22-Q11_pi/3} then $Q_{22}(\pi/3,t)-Q_{11}(\pi/3,t)>0.741+0.1-0.84=0.001>0$, which shows that $Q_{22}(\pi/3,t)>Q_{11}(\pi/3,t)$. This completes the proof of~\eqref{item:5}.
\end{proof}
\let\qed\relax
\end{proof}

\section*{Declaration of interests}

The authors do not work for, advise, own shares in, or receive funds from any organization that could benefit from this article, and have declared no affiliations other than their research organizations.

\printbibliography

\end{document}