\section{Introduction}
The sensing of individual electron charges confined to quantum dots (QD) forms the basis of many important quantum mechanical measurements in solid state systems such as spin readout, transfer and coherent manipulation~\cite{dicarlo2004,braakman2013,cao2012,watson2015}. As the complexity of coupled QD devices increases so does the electronic tunnelling processes that occur within them~\cite{schleser2004,PhysRevLett.112.176803,watson2015b}. In particular, for multi-QD systems the tunnelling of electrons from dot-to-dot will often occur in multiple stages resulting in non-trivial tunnelling paths. Establishing the path an individual electron takes is important when quantum information is being stored on specific particles within a multi-QD system~\cite{morello2010,buch2013,A.:2016fk}.

In any multi-QD architecture, measurement of the steady state charge occupation via a charge stability map~\cite{PhysRevB.67.161308,petta2004manipulation,PhysRevLett.94.196802,A.:2016fk,taylor2007} is a prerequisite for more complex measurements such as spin readout~\cite{elzerman2004,morello2010,buch2013}. However, these maps only present the long-term charge equilibrium of the system. This is because a charge stability measurement reveals information about the chemical potential of quantum dots with respect to their electron reservoir, not with respect to one another leaving the tunnelling path taken by an individual electron through the system unknown, see Fig.~\ref{fig:device}~\cite{hanson2007}. In order to determine the actual tunnelling path a particular electron has taken it is necessary to perform time resolved charge sensing~\cite{kung2009} thereby measuring the non-equilibrium dynamics of the QD system~\cite{house2011}.

The smallest system where non-trivial electron tunnelling paths can occur is a triple quantum dot (TQD). Since their first realisation, a host of new physics has been investigated using the TQD~\cite{PhysRevLett.75.705} including new tunnelling regimes~\cite{watson2014,PhysRevB.73.235310,PhysRevLett.112.176803}, novel spin-blockade effects~\cite{busl2013}, and complex coherent spin physics~\cite{gaudreau2012,laird2010,A.:2016fk} for use in quantum information. In this paper we investigate the non-equilibrium behaviour of single-electron tunnelling in a TQD fabricated using precision placed donor atoms in silicon~\cite{fuechsle2012,watson2014,gorman2015b}. We use real-time charge sensing to determine the preferred tunnelling paths of electrons over a range of electrostatic potentials. Importantly we develop a Anderson-Hubbard model that describes the various electron pathways through multi-QD structures.

The paper is set out as follows: In section~\ref{sec:the} we outline the model used to describe the electronic tunnelling pathways in multi-QD systems. Importantly, the device has been designed so that dot-dot tunnel rates are much greater than dot-reservoir tunnel rates. This means that the individual dot-reservoir rates are the limiting inputs for this model. From it we show we can accurately determine the tunnelling pathways under any electrostatic configuration. In section~\ref{sec:exp} we describe the device design and fabrication used to demonstrate our real-time charge sensing technique, the results of which are presented in section~\ref{sec:res}. Finally in section~\ref{sec:dis} we discuss our results in the context of an increasingly complex QD architecture and show how our results can be used to elucidate tunnelling pathways in multi-electron systems.

\section{Theory of Multi-Dot Tunnelling}
\label{sec:the}

The most natural way to describe the tunnelling dynamics of electrons from multiple QDs is to use a network of coupled `sites', where each site has both an intra- and inter-site energy $\epsilon$ and $U$ respectively. For $M$ QDs each hosting at most $n{=}1$ electrons coupled to a single reservoir the system is embodied by a modified Hubbard Hamiltonian given by,
\begin{equation}
H = \sum^{M}_j (\epsilon_j - \mu_j) \hat{n}_j + \sum_{j{\neq}k} \frac{U_{j,k}}{2} \hat{n}_j \hat{n}_k + \sum_{i} \epsilon_i\hat{n}_i,
\end{equation}
where $\epsilon_j$ is the detuning, $\mu_j$ is the chemical potential of QD $j$, the number operator for QD $j$ is $\hat{n}_j$ and $U_{j,k}$ is the inter-Coulomb repulsion between QDs $j$ and $k$. The last term represents the electronic states of the reservoir. In this work we consider QDs on a linear graph with nearest-neighbour coupling, however, the model can be extended to an arbitrarily complex graph of QDs and reservoirs. The detunings $\epsilon_j$ are determined by voltages applied to gates, $V_g$, surrounding the device,
\begin{equation}
\label{eq:lever}
\epsilon_j = \sum_{g} \alpha_{j,g} V_g,
\end{equation}
where $\alpha_{j,g}$ is the conversion factor (or lever arm) from applied voltage to energy from gate, $g$ to QD $j$~\cite{vanderwiel2003}. The complete electrostatic configuration is therefore deduced by considering all QD detunings which we will label $\Lambda{=}\{\epsilon_1,... \epsilon_j,... \epsilon_M\}$. We omit the coherent coupling terms between the QDs and instead assume that tunnelling between them is completely incoherent. This assumption is valid over the energy scale investigated in this work ($\sim10$~meV) which is much larger than regions where coherent tunnelling can occur ($\sim10$~$\mu$eV) i.e. at inter-dot transitions.

\begin{figure}[t]
\begin{center}
\includegraphics[width=1\columnwidth]{Figure_1.png}
\end{center}
\caption{{\bf A triple quantum dot in Si:P with an adjacent single-electron-transistor used as a charge sensor.} a) A scanning tunnelling micrograph of three small Si:P QDs incoherently tunnel coupled to a larger QD which itself is coupled to source ($S$) and drain ($D$) leads and acts as a sensing single electron transistor (SET). A gate $G_\text{SET}$ is used to control the SET while gates $\{G_{1},G_{2},G_{3}\}$ are used to control the QDs $\{1,2,3\}$. Inset shows a close up image of the three QDs. From the lithographic area it is estimated that they contain $\sim5$ donors each. The distances between the QDs as well as their individual distances to the SET are shown in nanometers. b) The conductance of the SET as a function of the gates $\{G_\text{1},G_\text{3}$\} for $G_\text{2}=0.9$~V and $G_\text{SET}=0.5$~V. Breaks in the Coulomb blockade peaks (lines running at $-45^\circ$) of the SET show where a charge transition from one of the QDs to the SET occurs. The inset provides a guide for where the charge transitions occur. The equivalent charge numbers on the QDs are given by $(n_1, n_2, n_3)$. c) A schematic representation of the modified Hubbard model showing the energies (not to scale) of different charge configurations at a detuning given by the star, also shown in (b). The incoherent tunnel rate $\Gamma_{s,s'}$ couples two charge configurations $s$ and $s'$. }
\label{fig:device}
\end{figure}

Finally, to incorporate the incoherent tunnelling in the system we transform the system from Hilbert space to Liouville space~\cite{gorman2015a} and introduce the incoherent rates $\Gamma_{s,s'}$ for the charge transition $s{\rightarrow}s'$. By doing this we can neglect the many charge states of the reservoir and consider only the $2^M$ possible charge states over the $M$ QDs. Importantly, the rates $\Gamma_{s,s'}$ represent direct transitions only, i.e. dot-to-dot or dot-to-reservoir, and do not represent any indirect tunnelling of an electron involving multiple tunnelling events. A schematic representation of the model is shown in Fig.~\ref{fig:device}c.

These incoherent tunnel rates $\Gamma_{s,s'}$ follow a Fermi-distribution dependent on the difference in energy $\Delta E_{s,s'}$ between the charge states $s$ and $s'$,
\begin{equation}
\Gamma_{s,s'}(\Lambda){=}\begin{cases}
\frac{\eta_{j,k}}{1 + \textnormal{exp}[\Delta E_{s,s'}(\Lambda)/k_B T]}, & \Delta N = 0 \\ \\
\frac{\gamma_{j} - \beta(\Delta E_{s,s'}(\Lambda)-\mu_{\text{\tiny{RES}}})}{1 + \textnormal{exp}[(\Delta E_{s,s'}(\Lambda)-\mu_{\text{\tiny{RES}}})/k_B T]}, & \Delta N = 1,
\end{cases}
\label{eq:rate}
\end{equation}
\noindent where $N$ is the total number of electrons across all QDs; $\eta_{j,k}$ is the tunnel rate from QD $j{\rightarrow}k$ for a change $\Delta N{=}0$ and $\gamma_j$ is the tunnel rate from QD $j$ to the reservoir for $\Delta N{=}1$~\cite{gorman2015b}. Here we denote the chemical potential of the reservoir as $\mu_{\text{\tiny{RES}}}$. The thermal energy is $k_B T$ at a temperature, $T$, and the energy difference between two charge states, $s$ and $s'$, is given by, $\Delta E_{s,s'}(\Lambda){=}E_s(\Lambda){-}E_{s'}(\Lambda)$. To account for an observed linear increase in tunnel rates from the QDs to the reservoir as a function of detuning we include a phenomenological dimensionless factor, $\beta$, that increases the tunnel rate for large values of $\Delta E_{s,s'}(\Lambda){-}\mu_{\text{\tiny{RES}}}$. From a cut in the experimental data shown in Fig.~\ref{fig:theoryvsexpt}e as the white dashed line, we estimate this value to be, $\beta{=}2{\times}10^{-3}$. Following Fermi's golden rule, this factor accounts for the change in matrix element coupling the QD and the reservoir as the QD is detuned from the $\mu_j{=}\mu_{\text{\tiny{RES}}}$ condition~\cite{schleser2004}.

Using this model the tunnel rate from an initial excited charge state $\rho_i$ to the ground state $\rho_{gs}$ for the electrostatic configuration $\Lambda$ is determined by solving the Lindblad master equation~\cite{contreras-pulido2014,lindblad1976},
\begin{equation}
\frac{d\rho}{dt} = (\mathcal{L}_c + \mathcal{L}_{\Gamma})\rho,
\end{equation}
where $\mathcal{L}_c = i(\mathbb{I} \otimes H - H \otimes \mathbb{I})$ is the coherent time evolution and the incoherent term is given by,
\begin{multline}
\mathcal{L}_{\Gamma} = \sum_{s,s'} \frac{\Gamma_{s,s'}(\Lambda)}{2} (2 L_{s,s'} \otimes L_{s,s'} - L_{s,s'}' L_{s,s'} \otimes \mathbb{I}\\ - \mathbb{I} \otimes L_{s,s'}' L_{s,s'}),
\end{multline}
where $L_{s,s'} = \ketbra{s}{s'}$. The thermal ground state of the system at $\Lambda$ is then determined by,
\begin{equation}
\rho_{gs}(\Lambda) = \frac{1}{Z} e^{-H(\Lambda)/k_B T},
\end{equation}
\noindent where, $\rho_{gs}(\Lambda)$ is the density operator for the thermal ground state and $Z{=}\textnormal{Tr}(e^{-H(\Lambda)/k_B T})$ is the partition function of the system. For our model we assume an electron temperature of $T{=}200$~mK, which is a typical value for our devices measured in a dilution fridge~\cite{buch2013,weber2014,gorman2015b}. After allowing the system to evolve from $\rho_i{\rightarrow}\rho_{gs}$, the tunnel time, $1/\Gamma(\Lambda)$ is taken to be the time when the probability of the charge ground state $\rho_{gs}$ reaches $1-e^{-1}$. By repeating this procedure under different electrostatic configurations, $\Lambda$, we can produce a so-called \textit{tunnel rate map} for the multi-QD system. Unlike the canonical charge stability map, the tunnel rate map can reveal information about the position of QD chemical potentials with respect to one another over the available gate range.

\begin{figure*}
\includegraphics[width=0.75\textwidth]{Figure_2.png}
\caption{{\bf Comparison of theoretical and experimental tunnel rate maps.} a) The theoretically predicted tunnel rate map for the TQD device scaled by the highest dot-reservoir tunnel rate, $\gamma_2$. Solid magenta lines show the charge transitions visible in standard charge stability maps (see Fig.~\ref{fig:device}b) while the dashed and dotted black lines correspond to transitions that can only be obtained from the tunnel rate map. The inset shows a zoom-in around the (0,1,0) charge region. b) The points (i), (ii), (iii), and (iv) corresponding to those shown in (a) are examined in terms of the chemical potentials of the QDs and the SET Fermi level (shown by the grey shaded region on the right of each panel). The levels shown in red and green correspond to the chemical potentials of those charge states at the points shown by the red and green circles in (a). The dashed lines in (i) and (ii) represent the SET Fermi level and corresponds to the dashed equivalent line shown in (a). Likewise, the dotted lines in (iii) and (iv) represent the relevant chemical potentials of charge states also shown in (a). c) The experimental procedure used to produce a tunnel rate map involves a two-level pulse from the (1,1,1) charge region (red square) to the detuning position $\Lambda$ (green square). The position marked by the red square is stepped across the $V_{G1}{-}0.6V_{G3}$ direction, while the final detuning position shown by the green grid is stepped across the charge transitions (magenta lines). d) The SET current is monitored in time to detect any electron tunnelling during the two level pulse. All data (blue circles) is comprised of 200 repetitions of the two level pulse (black line). A measure of the tunnel rate $\Gamma(\Lambda)$ is obtained from an exponential fit (green line) during the first phase of the two level pulse. The second phase of the pulse indicates the reloading of the electrons into the (1,1,1) charge region. e) A tunnel rate is obtained for different final detuning positions (green marker in (c)) $\Lambda$. There is a good qualitative agreement between the theoretical tunnel rate map shown in (a) which predicts all of the features obtained from the measurement. The areas enclosed by the white lines indicate the regions in gate space where an electron initially on QD-$1$ can only tunnel directly to the SET. The dashed white line indicates the cut used to calculate $\beta$ for Eq.~\ref{eq:rate}.}
\label{fig:theoryvsexpt}
\end{figure*}

\section{Experiment}
\label{sec:exp}
The device studied in this work is a TQD fabricated using scanning tunnelling microscopy hydrogen lithography. The methods of this fabrication technique have been reported in detail previously~\cite{fuechsle2007,simmons2008,fuhrer2009}. The TQD device shown in Fig.~\ref{fig:device} is comprised of three small QDs, labelled $1$, $2$, and $3$ (left, middle, and right respectively), consisting of $\sim$\,5~P donors each, determined by examining the extent of the exposed lithographic area~\cite{buch2013,weber2014}. These small donor clusters are tunnel coupled to a large quantum dot made up of $\sim 1000$P atoms and placed at a distance of 17--21~nm from it, see Fig.~\ref{fig:device}a. In turn this larger dot is coupled to source ($S$) and drain ($D$) leads allowing electrons to flow across its tunnel junctions on either side, therefore acting as a single-electron-transistor (SET) charge sensor~\cite{kastner1992,fuechsle2012}. The SET  has a charging energy of $5{\pm}1$~meV and is operated with a source-drain bias of 0.3~mV. The electrostatics of the SET are controlled by the gate $G_{SET}$, whereas the QDs are predominantly tuned using the gates $G_1$, $G_2$ and $G_3$. The SET island serves as the electron reservoir in this system and the incoherent coupling rates of the QDs $\{1,2,3\}$ are given by \{$\gamma_{1},\gamma_{2},\gamma_{3}$\}, respectively.

For the experiments presented in this letter $G_2$ is used as a global gate to shift the potential of all the QDs, while $G_1$ and $G_3$ are used to detune the potential of the QDs with respect to the SET. As such, the relevant lever arms, $\alpha_{j,g}$ for Eq.~\ref{eq:lever} are those of $G_1$ and $G_3$. A charge stability diagram of the TQD system, showing the SET conductance as a function of $G_1$ and $G_3$, is presented in Fig.~\ref{fig:device}b. Lines running at ${\sim}45^{\circ}$ in the data show the Coulomb blockade of the SET and breaks in these lines correspond to charge transitions between the SET and the three QDs. Due to the different capacitive coupling of the gates $\{G_1,G_3\}$ to each of the QDs, three distinct lines of SET breaks with different slopes are visible in this gate space. In addition a characteristic pentagon structure associated with the quadruple point of a TQD around (0,1,0) can be seen (see inset of Fig.~\ref{fig:device})~\cite{busl2013,braakman2013,watson2014}. We note that the absolute electron number are not known for this device; however, for the purpose of this work we assign the charge states shown in Fig.~\ref{fig:device}b where ($n_1,n_2,n_3$) represent the electron numbers on QDs $1$, $2$ and $3$, respectively. This does not affect the physics discussed in this work.

\section{Results}
\label{sec:res}
Using the model described in section~\ref{sec:the} and initialising the system in the $s{=}(1,1,1)$ charge configuration we obtain a theoretical tunnel rate map for this device shown in Fig.~\ref{fig:theoryvsexpt}a. The electrostatic configuration $\Lambda$ in gate space $V_{G1}{-}0.6V_{G3}$ vs $V_{G1}{+}V_{G3}$ was determined by calculating the lever arms, $\alpha_{j,g}$ using a capacitance modelling program~\cite{weber2012a}. This theoretical tunnel rate map  takes as an input the individual tunnel rates of the three QDs to the SET reservoir which, due to both their different distances from the SET and donor numbers, are given approximately by $\gamma_{1}{=}150$~Hz, $\gamma_2{=}1200$~Hz and $\gamma_3{=}1000$~Hz (obtained from experiment). In another work on the same device~\cite{gorman2015b} the tunnel couplings between dots 1 and 2 and between dots 2 and 3 were measured to be $5.5$~GHz and $2.2$~GHz respectively, which are much greater than any dot-reservoir rate. Thus, as an input for the model it is sufficient to use a incoherent dot-dot rate that is much larger than any dot-reservoir rate in order to obtain the correct system dynamics, here we assume $\eta_{1,2}{=}\eta_{2,3}{=}1000\gamma_{1}$. This argument is valid for any device with inter-dot separations of $\sim$~15nm in donor based systems, which typically results in tunnel coupling rates of the order $10^6$-$10^9$~Hz~\cite{weber2014,House2015ng}.

A tunnel rate map can be obtained experimentally following the procedure shown schematically in Fig.~\ref{fig:theoryvsexpt}c. At $t{=}0$ the device is initialised in the equivalent $(1,1,1)$ charge region, at the position labelled by the solid red square, from here an initial pulse is applied simultaneously to gates $G_1$ and $G_3$ to a position marked by the solid green square. During the first pulse phase, in cases where the ground state of the system contains a total electron number $<$3, one or more electrons will tunnel off the QDs to the SET reservoir during this pulse. Importantly, the length of this pulse, $t_1{=}30$~ms, is longer than any dot-to-reservoir tunnelling time, thus always allowing the system to reach its ground state charge configuration. A subsequent pulse is applied at time $t{=}t_1$ of length $t_2{=}30$~ms to reinitialise the charge state of the QDs into the ${(1,1,1)}$ charge configuration. Again, $t_2$ is much longer than any of the reservoir-to-dot tunnelling time. The same two level pulse sequence is repeated $200$ times and the average current of the SET is recorded over the duty cycle of one complete pulse sequence from $t{=}0$ to $(t_1{+}t_2)$, see Fig.~\ref{fig:theoryvsexpt}d.

In cases where at least one electron has tunnelled to the SET during the first phase of the pulse cycle the SET current will show an exponential change from which the experimental tunnel rate $\Gamma_{E}(\Lambda)$ is extracted. Note that, although theoretically multiple exponents of decay should be observed, in practice tunnelling will be dominated by the slowest rate, we therefore fit to a single exponential decay. The position of the first pulse is stepped across the charge stability region in segments shown by the subsequent green squares in Fig~\ref{fig:theoryvsexpt}c. In doing so, a tunnel rate can be deduced for every detuning position, $\Lambda$. Our measured tunnel rate map shown in Fig~\ref{fig:theoryvsexpt}e agrees very well with the theoretical map in Fig.~\ref{fig:theoryvsexpt}a. As well as showing the direct dot-reservoir charge transitions apparent in the standard charge stability map in Fig.~\ref{fig:device}b, the tunnel rate map reproduces features that arise from indirect tunnelling pathways.

As we show below, because the tunnel rates of the individual QDs to the SET reservoir have the relationship $\gamma_1{\ll}\gamma_3{<}\gamma_2$, in some circumstances it will be favourable for electrons to tunnel via other QDs as a means to reach the ground state in the shortest time possible, rather than tunnel directly to the SET.  Importantly, whether or not these types of tunnelling paths are allowed depends on the relative positions of the QD chemical potentials, which can be calculated using the constant interaction model~\cite{vanderwiel2003} (see Fig.~\ref{fig:theoryvsexpt}b). The positions of these chemical potentials can therefore be deduced directly from the tunnel rate map and we now discuss in detail the four points labelled in Fig.~\ref{fig:theoryvsexpt}a as (i-iv) that result in non-trivial tunnelling pathways.

At position (i) in Fig.~\ref{fig:theoryvsexpt}a the ground-state charge configuration is (0,1,1) and the chemical potentials of the QDs and SET at the position shown by the red circle are related by $\mu_1\text{\small{(1,1,1)}}{>}\mu_2\text{\small{(1,1,1)}}{>}\mu_{\text{\tiny{RES}}}{>}\mu_3\text{\small{(1,1,1)}}$ (where the subscript refers to an electron transition from this dot), see Fig.~\ref{fig:theoryvsexpt}b. Since the tunnel rates for QD-$1$ and QD-$2$ are related by $\gamma_1{<<}\gamma_2$ an electron is much more likely to tunnel to the SET from QD-$2$ first before the tunnelling of an electron from QD-$1$. The charge state will now be $(1,0,1)$, leaving the electron on QD-$1$ free to tunnel across to QD-$2$ because $\mu_1\text{\small{(1,0,1)}}{>}\mu_2\text{\small{(0,1,1)}}$. That is, the system will effectively make the following state transitions: $(1,1,1){\rightarrow}(1,0,1){\rightarrow}(0,1,1)$. The reduction in tunnel rates at the point indicated by the green circle indicates where $\mu_2\text{\small{(1,1,1)}} < \mu_{\text{\tiny{RES}}}$ where a slower tunnel rate results because an electron can only tunnel to the SET from QD-$1$ at this point. The boundary between these two regions (black dashed line) indicates where $\mu_2\text{\small{(1,1,1)}}{=}\mu_{\text{\tiny{RES}}}$.

The situation at position (ii) is almost identical to that at position (i), only here the ground-state charge configuration is (1,1,0) i.e. QD-$3$ takes the place of QD-$1$ at position (i). Again, the interface between the regions at (ii) $\mu_2\text{\small{(1,1,1)}}{=}\mu_{\text{\tiny{RES}}}$ is given by the black dashed line. At the point labelled with the red circle the charge transitions $(1,1,1){\rightarrow}(1,0,1){\rightarrow}(1,1,0)$ are made, whereas at the green circle the electron initially on QD-$3$ tunnels directly to the SET. It is worth noting that the same $\mu_2\text{\small{(1,1,1)}}{=}\mu_{\text{\tiny{RES}}}$ line can be mapped all the way through from the (1,1,0) to the (0,1,1) charge region.

At position (iii) the ground-state charge configuration is (0,1,0) and tunnelling from this region involves multiple tunnelling events over numerous possible pathways. Unlike at positions (i) and (ii) however, the tunnelling via different pathways depends on the relative \textit{inter-dot} chemical potentials. To illustrate this we will discuss one such tunnelling path in detail, the relevant energy diagram for which is shown in Fig.~\ref{fig:theoryvsexpt}b (iii). Here, $\mu_1\text{\small{(1,1,1)}}{>}\mu_2\text{\small{(1,1,1)}}{>}\mu_3\text{\small{(1,1,1)}}{>}\mu_{\text{\tiny{RES}}}$, i.e. any electron on any dot can in principle tunnel to the SET. In the case where QD-$3$ tunnels first, we are left with (1,1,0), and at the point shown by the red circle we have $\mu_2\text{\small{(1,1,0)}}{>}\mu_3\text{\small{(1,0,1)}}$ such that an electron will tunnel from QD-$2$ to QD-$3$. Finally, because $\mu_1\text{\small{(1,0,1)}}{>}\mu_2\text{\small{(0,1,1)}}$ the electron originally on QD-$1$ tunnels across to QD-$2$. The total electron movement in the system is: $(1,1,1){\rightarrow}(1,1,0){\rightarrow}(1,0,1){\rightarrow}(0,1,1){\rightarrow}(0,1,0)$. At the detuning position marked by the green circle however, we have $\mu_3\text{\small{(1,0,1)}}{>}\mu_2\text{\small{(1,1,0)}}$, which restricts the electron on QD-$2$ tunnelling over to QD-$3$ therefore forcing the electron on QD-$1$ to tunnel from there to the SET reservoir. This restriction reduces the observed tunnel rate because some of the tunnelling pathways will be limited by the slowest rate $\gamma_1$. A line separating two regions at (iii) indicates an alignment of the QD chemical potentials $\mu_2\text{\small{(1,1,0)}}{=}\mu_3\text{\small{(1,0,1)}}$.

The last position we discuss is far outside of the TQD quadruple point itself and is shown by position (iv) in Fig.~\ref{fig:theoryvsexpt}a. Here the ground state is (0,0,0), and between the green and red circles a transition can be seen where the two QD chemical potentials $\mu_1\text{\small{(1,0,0)}}{=}\mu_2\text{\small{(0,1,0)}}$, i.e. where an inter-dot electron transition can occur. Here, all three electrons can tunnel off to the SET, but the order in which they do so depends on the relative position of the QD chemical potentials with respect one another. At the position labeled by the green circle where $\mu_1\text{\small{(1,0,0)}}{<}\mu_2\text{\small{(0,1,0)}}$ the electron on QD-$2$ and $3$ tunnel to the SET first, followed by the slow tunnelling of the QD-$1$ electron at the rate $\gamma_1$. However, at the red circle after the electron from QD-$2$ has tunnelled the electron now residing on QD-$1$ can tunnel across to QD-$2$ where it will finally tunnel to the SET but now at the much faster rate of $\gamma_2$. A similar effect cannot be seen at the equivalent transition between QD-2 and QD-3 where $\mu_3\text{\small{(0,0,1)}}{=}\mu_2\text{\small{(0,1,0)}}$ because the tunnel rates, $\gamma_2$ and $\gamma_3$, are not different enough.

It is worth noting that the same information can be obtained from a tunnel rate map for an arbitrary set of tunnel rates, $\gamma_i$ despite different preferred tunnelling pathways. However, the visibility of the transitions will be reduced as the difference between the tunnel rates approaches zero.

\begin{figure}
\begin{center}
\includegraphics[width=1\columnwidth]{Figure_3.png}
\end{center}
\caption{{\bf Time-resolved charge stability maps for a triple quantum dot.} Maps showing the conductance of the SET charge sensor at different times from $t{=}0.32$~ms to $t{=}16.9$~ms. The solid magenta lines show those transitions that are apparent from a standard stability map, whereas the dashed lines are only visible using time resolved charge sensing. a) The conductance through the SET charge sensor at 0.32~ms before any electrons have tunnelled off the QDs. There are no SET breaks corresponding to QD transitions. b) At 2.43~ms, both the electrons from QD-$2$ and $3$ have tunnelled to the SET creating breaks in the Coulomb blockade where $\mu_2$ and $\mu_3$ are equal to $\mu_{\text{\tiny{RES}}}$. Point (i) shows where $\mu_2\text{\small{(1,1,1)}}{=}\mu_{\text{\tiny{RES}}}$; whereas point (ii) indicates the $\mu_1\text{\small{(1,0,0)}}{=}\mu_2\text{\small{(0,1,0)}}$ condition. c) At 7.49~ms, the electron on QD-$1$ begins to tunnel to the SET, thus two transitions in the upper right corner of the map are visible at points (iii) and (iv). At point (iii) the electron is tunnelling straight from QD-$1$ to the SET. The electron on QD-$1$ is also able to tunnel from the (1,0,0) charge state to the (0,0,0) on the bottom left of the map at point (iv). d) Finally, after a time 16.9~ms the map is equivalent to the standard charge stability map from Fig.~\ref{fig:device}b, as it shows all three of the QD transitions at equilibrium.}
\label{fig:gg_maps}
\end{figure}

To further reveal the non-equilibrium dynamics of electron movement across the TQD, we examine the instantaneous conductance through the SET charge sensor during the first pulse phase, between $t{=}0$ and $t_1$. These \textit{time resolved} charge stability maps are shown in Fig.~\ref{fig:gg_maps}. Figure~\ref{fig:gg_maps}a shows the instantaneous current through the SET at time $t{=}0.32$~ms. At this time no clear breaks in the SET Coulomb blockade peaks are apparent indicating that no electrons have tunnelled from the TQD system to the SET reservoir before this time, and therefore the system remains in the (1,1,1) charge configuration. In contrast, at $t{=}2.43$~ms, shown in panel (b), multiple SET breaks corresponding to charge transitions from the TQD to the SET reservoir can be seen. Not surprisingly, the line breaks observed at this time are associated with electron transitions from QD-$2$ and $3$ because they have the fastest tunnel rates.

Importantly, the map in Fig.~\ref{fig:gg_maps}b contains features seen in the standard charge stability map as well as some additional ones not seen in Fig.~\ref{fig:device}b that arise due to non-trivial tunnelling pathways. In particular an SET break corresponding to where the chemical potentials $\mu_2\text{\small{(1,1,1)}}{=}\mu_{\text{\tiny{RES}}}$ at point (i) in this map is clearly visible. This occurs since the electron on QD-$1$ has not had enough time to tunnel directly to the SET, yet the electron on QD-$2$ has tunnelled to the SET leaving space for the electron on QD-$1$ to tunnel across to this site (see discussion of (i) of Fig.~\ref{fig:theoryvsexpt}a above). Break (ii) in this map corresponds to where $\mu_1\text{\small{(1,0,0)}}{=}\mu_2\text{\small{(0,1,0)}}$ (see discussion of point (iv) from Fig.~\ref{fig:theoryvsexpt}a above).

The map shown in Fig.~\ref{fig:gg_maps}c is taken at a time where the electron on QD-$1$ has started to tunnel to the SET reservoir, and as a result, in addition to break (i) we see another transition in the top right corner of this map, break (iii), where $\mu_1\text{\small{(1,1,1)}}{=}\mu_{\text{\tiny{RES}}}$.  Since we see both transitions at this time the system must have a non-zero probability of tunnelling via two pathways near this detuning region. In addition, another transition line corresponding to tunnelling from QD-$1$ is observed, at point (iv). Finally, after all the electrons have tunnelled at 16.9~ms in panel (d) the stability map of the TQD is fully recovered.

\section{Discussion}
\label{sec:dis}
In this work we consider the complex electron tunnelling pathways to a reservoir that arise in coupled QD systems. We have presented a detailed theoretical model that captures all of the non-equilibrium charge dynamics, and one that predicts very well the observed experimental signatures of non-trivial tunnelling processes. Although we consider only three QDs with nearest neighbour couplings, the model and results we present are applicable to any size and form of graph, including multiple QDs coupled to multiple reservoirs~\cite{seo2013,A.:2016fk,thalineau2012}. In these cases, one must add additional tunnel rates for individual reservoirs to Eq.~\ref{eq:rate}.  We also note that a reversal of our protocol, i.e. the \textit{loading} of electrons from a reservoir to the dots, will reveal an equivalent map of interdot chemical potentials.

Our experimental method relies upon individual QDs having different tunnel rates to distinguish between different tunnelling pathways. For donor based architectures, this caveat is easily fulfilled due to the sensitive dependency of tunnel rates on dot-reservoir distances (sensitive to the order of the order of $1$~nm displacements). In gate defined QDs, incoherent tunnel rates can be tuned using external gate biases such that a sufficiently large difference between individual QDs may be attained~\cite{PhysRevLett.75.705}. However, any measurable difference in tunnel rate is sufficient in order to establish a tunnel rate map.

The time resolved charge stability maps presented in Fig.~\ref{fig:gg_maps} give a direct image of the inter-dot chemical potentials, which can later be used to elucidate the complex tunnelling pathways that occur in these systems. Importantly, electron pathways from dot-to-reservoir can be dependent on the inter-dot chemical potential structure, both deep inside a stable charge region as well as near dot-reservoir transitions. Understanding and controlling electron pathways is a vital component of multi-electron physics, and to this end tunnel rate maps provide an important characterisation tool for solid-state quantum information processing. For example, some spin readout schemes require the transfer of an electron from one QD to another before readout can be carried out~\cite{Veldhorst2015zt,A.:2016fk}. In these cases it becomes imperative to characterise the movement of specific electrons under all electrostatic conditions in order that the location and movement of quantum information in the system can tracked and processed correctly.

\begin{acknowledgments}
This research was conducted by the Australian Research Council Centre of Excellence for Quantum Computation and Communication Technology (project no. CE110001027) and the US National Security Agency and US Army Research Office (contract no. W911NF-08-1-0527). M.Y.S. acknowledges an Australian Research Council Laureate Fellowship.
\end{acknowledgments}


\section{Introduction}
The suppression of quarkonia (bound states of a heavy quark and its anti-quark) in ultra relativistic heavy ion collisions is one of the most prominent probes used to investigate and quantify the properties of the quark gluon plasma (QGP). The in-medium dissociation probability of the different quarkonium states could provide an estimate of the temperature of the system since the dissociation is expected to take place when the medium reaches or exceeds the critical temperature of the phase transition ($T_{\rm c}$), depending on the binding energy of the quarkonium state. For charmonium ($c\overline c$) states, the J/$\psi$ is likely to survive significantly above $T_{\rm c}$ (1.5 - 2 $T_{\rm c}$) whereas $\chi_{\rm c}$ and $\psi(2\rm S)$ melt near $T_{\rm c}$ (1.1 - 1.2 $T_{\rm c}$)~\cite{satz,satz2}. At LHC energies, due to the large increase of the $c\overline c$ production cross-section with the collision energy, there is a possibility of J/$\psi$ production via recombination of ${c}$ and $\overline {c}$. Thus, the observation of J/$\psi$ production in nucleus-nucleus collisions via recombination also constitutes an evidence of QGP formation. The study of the $\psi(2\rm S)$ production, due to its different binding energy, is complementary to that of the J/$\psi$ and it may also be useful for the evaluation of the temperature of the medium. The $\psi(2\rm S)$-to-J/$\psi$ cross-section ratio is predicted to be very sensitive to the details of the recombination mechanism. Experimentally this ratio is interesting as most of the systematic uncertainties cancel, with the remaining systematic uncertainties being only due to the signal extraction and the efficiency evaluation. The pp results for the charmonium provide a baseline for the nuclear modification factor of charmonium production in \mbox{p-Pb} and \mbox{Pb-Pb} collisions. The study of charmonia in \mbox{p-Pb} collisions can be used as a tool for a quantitative investigation of the cold nuclear matter (CNM) effects including various mechanisms such as gluon shadowing, $c\overline c$ break-up via interaction with nucleons, initial/final state energy loss, relevant in the context of studies of the strong interaction. The region of very small $x$ is accessible at the LHC and therefore strong shadowing and coherent energy loss effects are expected.

\section{ALICE detector and data samples}
The ALICE Collaboration has studied $\psi$(2S) production through its dimuon decay channel, in the Muon Spectrometer which covers the pseudorapidity range $-$4 $< \eta <$ $-$2.5. The ALICE detector is described in detail in~\cite{jinst}. The pp analysis has been performed on a triggered event sample corresponding to an integrated luminosity of $\mathcal{L}^{\rm pp}_{\rm int}$ = 1.35 $\pm$ 0.07 pb$^{-1}$ in the rapidity interval $2.5 < y_{\rm lab} < 4$ at $\sqrt{s}$ = 7 TeV. The p-Pb data have been collected at $\sqrt{s_{\rm NN}}$ = 5.02 TeV under two different configurations, inverting the direction of the p and Pb beams. In this way both forward rapidity $2.03 < y_{\rm cms} < 3.53$ ($\mathcal{L}^{\rm pPb}_{\rm int}$ = 5.01 $\pm$ 0.19 nb$^{-1}$) and backward rapidity $-4.46 < y_{\rm cms} < -2.96$ ($\mathcal{L}^{\rm Pbp}_{\rm int}$ = 5.81 $\pm$ 0.18 nb$^{-1}$) could be accessed, with the positive $y$ defined in the direction of the proton beam. Finally, the Pb-Pb analysis has been performed at $\sqrt{s_{\rm NN}}$ = 2.76 TeV ($\mathcal{L}^{\rm PbPb}_{\rm int}$ = 68.8 $\pm$ 0.9 $\mu$b$^{-1}$) in the rapidity region $2.5 < y_{\rm lab} < 4$.
\section{Results}

\subsection{pp collisions}
Fig.~\ref{fig1} shows the inclusive differential production cross-sections of $\psi$(2S) as a function of $p_{\rm T}$ and $y$~\cite{epjc74}. The result on $p_{\rm T}$ differential  cross-section is consistent with the LHCb measurement~\cite{3epjc72} in the same rapidity interval. This is the first measurement of $\psi$(2S) differential cross-sections as a function of $y$ in pp collisions at $\sqrt{s} = 7$ TeV.

\begin{figure}[ht]
\includegraphics[scale=0.38]{2014-Sep-08-PSI_PT.pdf}
\includegraphics[scale=0.38]{2014-Sep-08-PSI_Y.pdf}
\caption{\label{fig1}Inclusive differential production cross-sections of $\psi$(2S) as a function of $p_{\rm T}$ (left) and $y$ (right).}
\end{figure}

The inclusive $\psi$(2S) to J/$\psi$ cross-section ratio was measured as a function of $p_{\rm T}$ and $y$ as shown in Fig.~\ref{fig2}. A clear $p_{\rm T}$  dependence can be observed, consistent with the one measured by LHCb~\cite{3epjc72}. No strong $y$ dependence is visible in the $y$ range covered by the ALICE muon spectrometer.

\begin{figure}[ht]
\includegraphics[scale=0.38]{2014-Sep-08-RATIO_PT.pdf}
\includegraphics[scale=0.38]{2014-Sep-08-RATIO_Y.pdf}
\caption{\label{fig2}Inclusive $\psi$(2S) to J/$\psi$ cross-section ratio as a function of $p_{\rm T}$ (left) and $y$ (right).}
\end{figure}

\subsection{p-Pb collisions}
The production cross section of $\psi$(2S) in p-Pb is compared to the J/$\psi$ one and to the corresponding quantities in pp collisions at $\sqrt{s}$ = 7 TeV (no LHC results are available at $\sqrt{s}$ = 5.02 TeV) using the $\psi$(2S) to J/$\psi$ ratio and the double ratio [$\sigma_{\rm \psi(2S)}/\sigma_{\rm J/\psi}]_{\rm pPb}$/[$\sigma_{\rm \psi(2S)}/\sigma_{\rm J/\psi}]_{\rm pp}$~\cite{jhep12,arnaldi} as shown in Fig.~\ref{fig3}. The pp ratios are significantly higher than those for p-Pb and Pb-p. The double ratio is compared with the corresponding measurement by the PHENIX Collaboration at mid-rapidity at $\sqrt{s_{\rm NN}}$ = 0.2 TeV~\cite{ada13}. Within uncertainties, a similar relative $\psi$(2S) suppression is observed by the two experiments.

\begin{figure}[ht]
\includegraphics[scale=0.38]{Fig2_Psi2S_JPsi_ratio_integrated.pdf}
\includegraphics[scale=0.38]{Fig3_Psi2S_JPsi_doubleratio_integrated_PhenixBlack.pdf}
\caption{\label{fig3}Left: $\psi$(2S) to J/$\psi$ cross-section ratio compared to the pp results at $\sqrt{s}$ = 7 TeV. Right: the double ratio compared to the PHENIX result~\cite{ada13}.}
\end{figure}

Since no result on cross-section of $\psi$(2S) is available at $\sqrt{s}$ = 5.02 TeV in pp collisons, the nuclear modification factor of $\psi$(2S) is obtained by combining the J/$\psi$ $R_{\rm pPb}$~\cite{jhep1402} and the double ratio, as $R_{\rm pPb}^{\rm \psi(2S)} =$$ R_{\rm pPb}^{\rm J/\psi}$$\times$$(\sigma_{\rm pPb}^{\rm \psi(2S)}/\sigma_{\rm pPb}^{\rm J/\psi})$$\times$$(\rm \sigma_{\rm pp}^{\rm J/\psi}/\sigma_{\rm pp}^{\rm \psi(2S)})$, assuming that the ratio in pp collisions does not depend on $\sqrt{s}$~\cite{jhep12}. In Fig.~\ref{fig4}, $R_{\rm pPb}^{\psi(2S)}$ is compared with $R_{\rm pPb}^{J/\psi}$ and also with theoretical calculations based on nuclear shadowing~\cite{ijmp}, coherent energy loss or both~\cite{jhep1303}. The suppression of $\psi$(2S) production is stronger than that of J/$\psi$ and reaches a factor of 2 with respect to pp. Since the kinematic distributions of gluons producing the J/$\psi$ or the $\psi$(2S) are rather similar and since the coherent energy loss does not depend on the final quantum numbers of the resonances, the same theoretical calculations hold for both J/$\psi$ and $\psi$(2S). Theoretical models predict a $y$ dependence which are in reasonable agreement with the J/$\psi$ results but no model can describe the $\psi$(2S) data. These results show that other mechanisms must be invoked in order to describe the $\psi$(2S) suppression in p-Pb collisions.

\begin{figure}[ht]
\centering
\includegraphics[scale=0.38]{Fig4_Psi2S_RpA_integrated.pdf}
\caption{\label{fig4}$\psi$(2S) $R_{\rm pPb}$ versus $y$ compared to the J/$\psi$ $R_{\rm pPb}$ and theoretical models.}
\end{figure}

The $R_{\rm pPb}$ is also computed as a function of $p_{\rm T}$ both at backward and forward $y$ and the results are shown in Fig.~\ref{fig5}. At both rapidities, the $R_{\rm pPb}^{\rm \psi(2S)}$ shows a strong suppression with a slightly more evident $p_{\rm T}$ dependence at backward-$y$. The $\psi$(2S) is more suppressed than the J/$\psi$, as already observed for the $p_{\rm T}$-integrated result. Theoretical calculations are in fair agreement with the $R_{\rm pPb}^{\rm J/\psi}$ but clearly overestimate the $R_{\rm pPb}^{\rm \psi(2S)}$ behaviour. The calculations from the comover model~\cite{ferre} shows that the interaction with comovers, mostly at play in the backward region, is able to explain the stronger $\psi$(2S) suppression.

\begin{figure}[ht]
\includegraphics[scale=0.38]{Fig7a_p-Pb_Psi2S_RpA_differential.pdf}
\includegraphics[scale=0.38]{Fig7b_Pb-p_Psi2S_RpA_differential.pdf}
\caption{\label{fig5}$p_{\rm T}$ dependence of the $\psi$(2S) $R_{\rm pPb}$ compared to the J/$\psi$ $R_{\rm pPb}$ and theoretical calculations in the forward (left) and backward (right) rapidity region.}
\end{figure}

Finally, the $\psi$(2S) production is studied as a function of the collision event activity both at backward and forward $y$~\cite{arnaldi}, as shown in Fig.~\ref{fig6}. The event activity determination is described in details in~\cite{toia}. Since the centrality determination in p-Pb collisions can be biased by the choice of the estimator, the nuclear modification factor is, in this case, named $Q_{\rm pPb}$~\cite{toia}. The $\psi$(2S) $Q_{\rm pPb}$ shows a strong suppression, which increases with increasing event activity, and is rather similar in both the forward and the backward $y$ regions. The J/$\psi$ $Q_{\rm pPb}$ shows a similar decreasing trend at forward-$y$ as a function of the event activity. On the contrary, the J/$\psi$ and $\psi$(2S) $Q_{\rm pPb}$ patterns observed at backward-$y$ are rather different, with the $\psi$(2S) significantly more suppressed for large event activity classes.

\begin{figure}[ht]
\includegraphics[scale=0.38]{2014-May-14-QpA_JPsi_Psi2S_EvActivity.pdf}
\includegraphics[scale=0.38]{2014-May-14-QAp_JPsi_Psi2S_EvActivity.pdf}
\caption{\label{fig6}$\psi$(2S) $Q_{\rm pPb}$ versus event activity compared to the J/$\psi$ $Q_{\rm pPb}$ in the forward (left) and backward (right) rapidity region.}
\end{figure}

\subsection{Pb-Pb collisions}
In order to study the suppression of $\psi$(2S) relative to J/$\psi$ in Pb-Pb collisions, the double ratio [$\sigma_{\rm \psi(2S)}/\sigma_{\rm J/\psi}]_{\rm PbPb}$/[$\sigma_{\rm \psi(2S)}/\sigma_{\rm J/\psi}]_{\rm pp}$ has been measured as a function of centrality in two $p_{\rm T}$ intervals (0 $<$ $p_{\rm T}$ $<$ 3 GeV/$\it{c}$ and 3 $<$ $p_{\rm T}$ $<$ 8 GeV/$\it{c}$)~\cite{arnaldi13} and has been compared with the results from the CMS Collaboration~\cite{prl113} as shown in Fig.~\ref{fig7}. Limited $\psi$(2S) statistics does not allow any firm conclusion about the centrality dependence of this ratio and the comparison with CMS is not straightforward due to the different kinematic coverage.

\begin{figure}[ht]
\centering
\includegraphics[scale=0.40]{2015-Sep-22-DoubleRatio_v12_Derived_ComparisonWithCMS_v3.pdf}
\caption{\label{fig7}Double ratio [${\rm \psi(2S)}/{\rm J/\psi}]_{\rm PbPb}$/[${\rm \psi(2S)}/{\rm J/\psi}]_{\rm pp}$ as a function of centrality in two $p_{\rm T}$ intervals compared to CMS measurements~\cite{prl113}.}
\end{figure}

\section{Summary}

In summary, the ALICE Collaboration has studied the inclusive $\psi$(2S) production in pp, p-Pb and Pb-Pb collisions at the LHC. The $\psi$(2S) production cross-section and the $\psi$(2S) to J/$\psi$ cross-section ratio have been obtained as a function of $p_{\rm T}$ and $y$ in pp collisions. The $p_{\rm T}$ differential results are in good agreement with LHCb measurements. In p-Pb collisions the results show that $\psi$(2S) is significantly more suppressed than J/$\psi$ in both rapidity regions. This observation implies that initial state nuclear effects alone cannot account for the modification of the $\psi$(2S) yields, as also confirmed by the poor agreement of the nuclear modification factor of $\psi$(2S) with models based on shadowing and/or energy loss. Interaction with comovers is able to explain the $\psi$(2S) suppression. The final state interaction with the hadronic medium could also provide a possible explanation for the stronger $\psi$(2S) suppression~\cite{rapp}. Limited statistics prevent to make definitive conclusions on $\psi$(2S) production in Pb-Pb collisions.


\section{Introduction}

The original idea of the R--matrix theory was introduced by Kapur and Peierls \cite{Kapur}, to remove unsatisfactory reliance on perturbation theory in nuclear reactions. It was a few years later that
Wigner simplified the idea to the formulation of R--matrix theory in which all expressions are energy dependent \cite{Wigner-Physrev-70,Wigner-Physrev-70p,Wigner-phys72}.
Although the R--matrix theory was originally developed for treatment of nuclear resonances \cite{Wigner-phys72,Lane-Mod}, it can also be used to describe all types of reaction phenomena and can be considered as an elegant method to solve the Schr\"{o}dinger equation. These extensions especially became possible after the work of Bloch \cite{Bloch}, by introducing a singular operator between internal and external regions.

The solution of the Schr\"{o}dinger equation for the bound states of few--quarks in momentum space is numerically difficult to handle, because the confining part of the potential leads to singularity at small momenta. To overcome this problem we have successfully used a regularized form of the quark--antiquark \cite{Hadizadeh_AIP1296} and recently diquark--antidiquark ($D\bar D$) \cite{Hadizadeh-PLB753,Hadizadeh-submitted} interactions to solve the Lippmann--Schwinger (LS) equation in momentum space and calculate the mass spectrum of heavy quarkonia and tetraquarks.
To this aim one can keep the divergent part of the potential fixed after exceeding a certain distance, called the regularization cutoff. This procedure creates an artificial barrier and the influence of tunneling barrier is manifested by significant changes in the energy
eigenvalues at small distances.

Recently we have shown that the homogeneous LS equation can be formulated in configuration space to study the heavy tetraquarks as a bound state of $D\bar D$ system \cite{Hadizadeh-EPJC75}.
The variational methods are also used to study the few--quark bound states \cite{Park_NPA925, Kamimura_PRA38, Brink_PRD57, Brac_ZPC, Vijande_PRD76, Vijande_EPJA19}. In order to obtain the variational energy, one must minimize the lowest eigenvalue with respect to the variational parameters after diagonalizing the Hamiltonian matrix. The variational energy can be obtained by differentiating the lowest eigenvalue with respect to the variational parameters.
The successful application of the R--matrix theory to describe the resonance and scattering states resulting from the interaction of particles or systems of particles, which can be nucleons, nuclei,
electrons, atoms and molecules \cite{Burke_2011,Descouvemont_RPP73}, motivated us to implement it to the few--quark bound states.
In this paper we have shown that R--matrix theory is an effective and efficient method to
study the bound states of few--quarks by solving the Schr\"{o}dinger equation in different partial wave channels.
The successful implementation of the R--matrix method to heavy few--quark systems, paves the path to accurately predict the masses of few--quark bound states composed of light quarks.

\section{R--matrix method for few--quark bound states in configuration space}

In our study for diquark, triquark, tetraquark and pentaquark systems, we have used a two--body picture by considering them as the bound states of quark--antiquark ($q\bar{q}$), diquark--quark ($Dq$), diquark--antidiquark ($D\bar{D}$) and diquark--antitriquark ($D\bar{T}$) systems, respectively.
The nonrelativistic bound state of any of these two--body systems with the pair relative distance ${\mathbf r}$ in a partial wave representation can be described by the Schr\"{o}dinger equation.
For simplicity of the notation we use $AB$ for representation of these systems, where $A$ and $B$ stand for any of the subsystems.
Since the $AB$ interaction $V(r)$ is a central force, thus the wave function of few-quark bound states consists of some type of radial function times a spherical harmonic function, i.e. $\Psi_{nlm}({\mathbf r}) = R_{nl}(r) Y_{lm}(\theta,\phi)$.
The radial function can be obtained by solution of the differential equation
\begin{eqnarray} \label{eq.Schrodinger}
&& \left ( H_{l}-E \right)R_{nl}(r) = \cr
&& \hskip-.5cm \left[-\dfrac{\hbar ^2}{2\mu} \Biggl (\dfrac{d^2}{dr^2}+\dfrac{2}{r}\dfrac{d}{dr}-\dfrac{l(l+1)}{r^2} \Biggr)+V(r)-E\right]R_{nl}(r)=0, \cr
\end{eqnarray}
where $E=m_{AB}-m_A - m_{B}$ is $AB$ binding energy ($m_{AB}$ is the mass of $AB$ system composed of two subsystems $A$ and $B$ with masses $m_A$ and $m_B$). $\mu = \dfrac{m_{A} m_{B}}{m_A +m_{B}}$ is the reduced mass and $l$ is the orbital angular momentum of relative motion of $AB$ system.
In R--matrix calculations, the radial wave function is considered in two internal ($r<r_c$) and external ($r>r_c$) regions as $R_{int}$ and $R_{ext}$, correspondingly. The parameter $a$ is large enough to be sure that the mass spectrum is independent of it.

 The internal wave function $R_{int}$ is defined as combination of the basis functions $u_i(r)$
\begin{equation} \label{Rintdirect}
 R_{int}=\sum_{i=1}^{N} c_i \, u_i(r),
\end{equation}
where
\begin{equation}
 u_i(r)=r^l \, e^{\left(-\dfrac{r^2}{2\rho_i^2} \right)}.
\end{equation}
In our study, the basis states parameters $\rho_i$ are
\begin{equation}
 \rho_i=\rho_0\lambda^{i-1}, \quad ~\rho_0=1.2\dfrac{r_c}{\lambda^{N-1}},\quad \lambda=1.3,\quad N=30.
\end{equation}

Since in the external region the interaction between diquark and antidiquark is fixed, the external wave function can be considered as modified spherical Bessel function of the second kind
\begin{equation}
R_{ext}=C_l \, k_l(\kappa r),
\end{equation}
where $C_l$ is a constant parameter and $\kappa=\sqrt{2 \mu\, (-E+V_0)}$.
Continuity of the internal and external wave functions and derivatives implies that
\begin{eqnarray} \label{Continuity}
 R_{int}(r_c) &=& R_{ext}(r_c), \cr
 R'_{int}(r_c)&=&R'_{ext}(r_c).
\end{eqnarray}

The Hamiltonian $H_{l}$ is not Hermitian in the internal region. To
avoid this the Bloch operator is defined as
\begin{equation}
 L(b)=\dfrac{\hbar^2}{2\mu} \, \delta \left(r-r_c \right) \, \left(\dfrac{d}{dr}-\dfrac{b}{r} \right).
\end{equation}
The dimensionless parameter $b$ is an arbitrary real constant. The delta function indicates that the Bloch operator is a surface operator acting only on $r = r_c$. The operator $H_{l} + L(b)$ is Hermitian, and therefore has a discrete spectrum in the finite region.
Using Eq. (\ref{Continuity}), the Schr\"{o}dinger equation in the internal region can be approximated by
\begin{equation}\label{eq.sch-bloch}
 \biggl ( H_{l}+L(b)-E \biggr ) \, R_{int}(r)=L(b) \, R_{ext}(r).
\end{equation}
It means the logarithmic derivative of the wave function is continuous at $r=r_c$.
By multiplying Eq. (\ref{eq.sch-bloch}) with $u_{i'}$
and integrating in the internal region,
we obtain the following equation to determine the unknown coefficients $c_i$
\begin{eqnarray}\label{eq2}
 \sum_{i=1}^N \, C_{i'i}(E,b) \, c_i &=&
 \langle u_{i'}|L(b)|R_{ext}\rangle \cr
 &=&
 \dfrac{\hbar^2 r_c}{2\mu} \, u_{i'}(r_c) \, \biggl ( aR'_{ext}(r_c)-bR_{ext}(r_c) \biggr),
\end{eqnarray}
where
\begin{equation} \label{eq.C}
 C_{i'i}(E,b)=\left \langle u_{i'} \left | H_{l}+L(b)-E \right | u_i \right \rangle.
\end{equation}

Solving Eq. (\ref{eq2}) for $c_i$ and substituting them into Eq. (\ref{Continuity}), i.e $R_{int}(r_c)=\sum_{i=1}^N c_i\, u_i(r_c) = R_{ext}(r_c)$, leads to
\begin{equation}
 R_{int}(r_c)=R(E,b) \, \biggl ( a\,R'_{ext}(r_c)-b\,R_{ext}(r_c) \biggr),
\end{equation}
where $R(E, b)$ is the R--matrix given by
\begin{equation} \label{eq.Reb}
 R(E,b)=\dfrac{\hbar^2 r_c}{2\mu} \sum_{i,i'=1}^N u_{i}(r_c) \, C^{-1}_{ii'}(E,b) \, u_{i'}(r_c).
\end{equation}

The wave function in the internal region is then given by
\begin{equation}\label{eq.Rint}
 R_{int}(r)=\dfrac{\hbar^2 r_c}{2\mu \, R(E,b)} \, R_{ext}(r_c) \sum_{i,i'=1}^N u_{i}(r) \, C^{-1}_{ii'}(E,b) \, u_{i'}(r_c).
\end{equation}

By choosing $b = \dfrac{a \, R'_{ext}(r_c)}{R_{ext}(r_c)} = \dfrac{\kappa r_c \, k'_l(\kappa r_c)}{k_l(k r_c)}$, where $k'_l(\kappa r_c) = \dfrac{dk_l}{d(\kappa r_c)}$, the right hand side of Eq. (\ref{eq2}) will be zero and consequently leads to the following Schr\"{o}dinger-Bloch equation
\begin{equation}
 \sum_{i=1}^{N} \left \langle u_{i'} \left |H_{l}+L \biggl (b(E) \biggr)-E \right | u_i \right \rangle c_i=0.
\label{eq.SB}
\end{equation}
Since $b$ depends on $\kappa$ or the $AB$ binding energy $E$ that we want to calculate, we have written $b \equiv b(E)$.
The equation (\ref{eq.SB}) can be written schematically as eigenvalue equation
\begin{equation} \label{eq.eigenvalue}
{\cal A}\, c = \lambda \, {\cal B} \, c
\end{equation}
where the matrix elements of ${\cal A}$ and ${\cal B}$ matrices can be obtained as
\begin{eqnarray}
 {\cal A}_{ij} &=& \left \langle u_{i} \left |H_{l}+L \biggl (b(E) \biggr) \right | u_j \right \rangle , \cr
 {\cal B}_{ij} &=& \left \langle u_{i} \left | \right. u_i \right \rangle.
\label{eq.AB-matrix}
\end{eqnarray}
Since the ${\cal A}$ matrix is energy dependent, the solution of the eigenvalue equation (\ref{eq.eigenvalue}) can be started by an initial guess for the energy $E$ and the search in the binding energy can be stopped when $\left |\dfrac{\lambda - E}{E} \right| \le 10^{-10}$.
In order to solve the Eq. (\ref{eq.SB}), we have discretized the continuous variable $r$ with Gauss-Legendre points using a hyperbolic--linear mapping \cite{Hadizadeh-PLB753}. To this aim we have transferred $[0,\infty)$ domain to $[0,1]\cup[1,2]\cup[2,15]\, \text{GeV}^{-1}$ using $75, 75$ and $50$ nodes in each subinterval, respectively.
It indicates that the parameter $r_c$, which divides the configuration space into internal and external regions, is chosen to be $r_c=15\, \text{GeV}^{-1}$ and we have numerically verified that the calculated masses of tetraquarks are independent of the regularization cutoff $a$.

By having the $AB$ binding energy and eigenvector $c$ from the solution of eigenvalue Eq. (\ref{eq.eigenvalue}), we can calculate the $AB$ internal wave function by equations
(\ref{eq.C}), (\ref{eq.Reb}) and (\ref{eq.Rint}). Of course, one can calculate the internal wave function directly using Eq. (\ref{Rintdirect}).
Using $AB$ wave function, we can evaluate the expectation value $\evdel{r}$ for $AB$ pair distance as
\begin{eqnarray}
\evdel{r} &=& \int_0^{\infty} dr\, r^3 R^2(r) \cr
&=& \int_0^{r_c} dr\, r^3 R_{int}^2(r) +  \int_{r_c}^{\infty} dr\, r^3 R_{ext}^2(r),
\end{eqnarray}
where the $AB$ radial wave function is normalized to 1, i.e. $\int_0^{\infty} dr\, r^2 R^2(r) = 1$.

\section{Results and Discussion} \label{results}

For numerical solution of the integral equation (\ref{eq.Schrodinger}) for $q \bar q$, $Dq$ and $D \bar D$ we have used the spin-independent interaction
\begin{equation}
 \label{eq.V_234}
V (r) = V_{Coul} (r) + V_{conf} (r),
\end{equation}
with the linear confining
\begin{equation}
 \label{eq.V2}
V_{conf} (r) = a \,r+b,
\end{equation}
and the Coulomb-like one-gluon exchange potential
\begin{equation}
 \label{eq.V2}
V_{Coul} (r) = \gamma \, \dfrac{F_{A} (r) F_{B} (r)}{r}, \quad \gamma = \dfrac{-4}{3} \alpha_s.
\end{equation}
$F_A$ and $F_{B}$ are the form factors of the subsystems $A$ and $B$, correspondingly, and have the following functional form
\begin{equation}
 \label{eq.F}
F (r)= 1-e^{\alpha r -\beta r^2}.
\end{equation}
The parameters of this model are fixed from the analysis of heavy quarkonia masses and radiative decays \cite{Galkin-SJNP44,Galkin-SJNP51,Galkin-SJNP55}.

\subsection{Heavy quarkonia}

For this first test of application of R--matrix method, we have solved the integral equation (\ref{eq.Schrodinger}) to calculate the mass spectra of heavy quarkonia, mesons consisting heavy quark and antiquark. We have used the linear confining plus coulomb potential of Eq. (\ref{eq.V_234}) with form factor $F(r)=1$.
The parameters of potentials are $a=0.18$ GeV$^2$, $b=-0.29$ GeV with $\alpha_s=0.47$ for charmonium ($m_c=1.56$ GeV) and $\alpha_s=0.39$ for bottomonium ($m_b=4.93$ GeV).
As we have shown in Table \ref{Table_qbarq}, our numerical results for masses of charmonium and bottomonium, obtained by R--matrix method are in excellent agreement with solution of Lippmann--Schwinger integral equation in momentum \cite{Hadizadeh_AIP1296} and configuration \cite{Faustov-IJMPA15} spaces and also with the experimental data \cite{Barnett-PRD54}.

\begin{table}[hbt]
\caption{The mass spectra of charmonium $\psi(c\bar{c})$ and bottomonium $\Upsilon(b \bar b)$ for the linear confining plus coulomb potential of Eq. (\ref{eq.V_234}) with form factor $F(r)=1$. The masses are given in GeV. The numbers in parentheses are the expectation value of the relative distance between quark--antiquark pair $\evdel{r}$ in units of fm.
}
\centering
\begin{tabular}{ccccccccccccccc}
\toprule
\multicolumn{1}{c}{State} && \multicolumn{4}{c}{$\psi(c\bar{c})$} &&  \multicolumn{4}{c}{$\Upsilon(b \bar b)$} \\
 \cline{3-6}   \cline{8-11}
 && R--matrix &  LS \cite{Hadizadeh_AIP1296}  & Faustov et al. \cite{Faustov-IJMPA15} & Exp. \cite{Barnett-PRD54}
 && R--matrix &  LS \cite{Hadizadeh_AIP1296}  & Faustov et al. \cite{Faustov-IJMPA15} & Exp. \cite{Barnett-PRD54}  \\ \hline
$1s$ && 3.062 (0.349)  & 3.062 & 3.068  &    3.0675 &&  9.421 (0.184) &    9.425 & 9.447  &   9.4604 \\
$1p$ && 3.529 (0.599) &  3.529 & 3.526  &    3.525  &&   9.910 (0.368) &    9.909 & 9.900   &   9.900 \\
$2s$ && 3.696 (0.734)  &  3.696 & 3.697  &   3.663   && 10.005 (0.453) &    10.006 & 10.012 &   10.023 \\
$1d$ && 3.832 (0.795) &  3.832 & 3.829  &    3.770  &&   10.158 (0.511) &    10.158 & 10.155 &    \\
$2p$ && 3.997 (0.920)  &  3.997 & 3.993  &              &&   10.263 (0.594) &    10.263 & 10.260 &  10.260 \\
$3s$ && 4.144 (1.040) &  4.144 & 4.144  &   4.159   &&  10.349 (0.669) &    10.350 & 10.353 &   10.355 \\
$2d$ && 4.238 (1.081) &  4.237 & 4.234  &               &&  10.451 (0.711) &  10.450 & 10.448 &    \\
$3p$ && 4.387 (1.222) &  4.384 & 4.383  &               &&   10.547 (0.786) &    10.546 & 10.544 &   \\
\bottomrule
\end{tabular}
\label{Table_qbarq}
\end{table}

\subsection{Heavy baryons}

In the next step we have calculated the masses of the ground state heavy baryons consisting of two light $(u; d; s)$ and one heavy $(c; b)$
quarks in the heavy--quark--light--diquark approximation.
The used diquark mass and form factor parameters are given in Table \ref{Table.pot_parameters}.
As we have shown in Table \ref{Table.baryon}, our numerical results for different heavy baryons, calculated by nonrelativistic R--matrix method, are in good agreement with relativistic and spin-dependent results of EFG~\cite{Ebert_PRD72} and also with MLW results \cite{Mathur_PRD66} of lattice nonrelativistic QCD.
The relative difference between our and Ebert, Faustov, and Galkin (EFG) group results is less than $1.2~ (4.8) \, \%$ for bottom (charm) baryons. Clearly the difference comes from the relativistic effects and also spin terms of the potential that we have ignored in our calculations.

\begin{table}[hbt]
\caption{The mass $m$ and form factor parameters $\alpha$ and $\beta$ of light--light and heavy--light diquarks. $q$ stands for up and down quarks, and $S$ and $A$ denote the scalar and axial vector diquarks.}
\centering
\begin{tabular}{cccccccccccc}
\toprule
quark  & Diquark  &  $m$ &$\alpha$ &$\beta$       \\
content & type & (GeV) & (GeV) & (GeV$^2$)
\\ \hline
\multirow{2}{*}{$qq$} & $S$ & 0.710 & 1.09 & 0.185  \\ \cmidrule{2-5}
& $A$ & 0.909 & 1.185 & 0.365  \\
\hline
\multirow{2}{*}{$qs$} & $S$ & 0.948 & 1.23 & 0.225  \\ \cmidrule{2-5}
& $A$ & 1.069 & 1.15 & 0.325  \\
\hline
$ss$ & $A$ & 1.203 & 1.13 & 0.280  \\
\hline
\multirow{2}{*}{$cq$} & $S$ & 1973& 2.55 &0.63 \\ \cmidrule{2-5}
& $A$ & 2.036& 2.51 &0.45 \\
\hline
\multirow{2}{*}{$cs$}& $S$ & 2091& 2.15 & 1.05 \\  \cmidrule{2-5}
 & $A$ & 2.158&2.12& 0.99 \\
\hline
\multirow{2}{*}{$bq$}& $S$ &  5359 &6.10 & 0.55 \\ \cmidrule{2-5}
& $A$ & 5.381& 6.05 &0.35 \\
\hline
\multirow{2}{*}{$bs$}& $S$ & 5462 & 5.70 &0.35 \\ \cmidrule{2-5}
& $A$ & 5.482 & 5.65 &0.27 \\
\bottomrule
\end{tabular}
\label{Table.pot_parameters}
\end{table}

\begin{table}[hbt]
\caption{Masses of the ground states of heavy charm and bottom baryons in units of GeV, calculated by nonrelativistic R--matrix method compared to results of EFG (from solution of relativistic LS equation), to MLW results (obtained from lattice nonrelativistic QCD) and also to the experimental date. The numbers in parentheses are the expectation value of the relative distance between light diquark and heavy quark pair $\evdel{r}$ in units of fm. The symbols $[]$ and $\{\}$ denote the scalar and axial vector diquarks.}
\centering
\begin{tabular}{cccccccccccc}
\toprule
Baryon & content & $I(J^P)$ & R--matrix & EFG~\cite{Ebert_PRD72} & MLW~\cite{Mathur_PRD66} & EXP PDG~\cite{Eidelman_PLB592}  \\ \hline
$\Lambda_c$ & $\{ss\}c$  & $0(\frac12^+)$ & 2.396 (0.49) &  2.297 & 2.290  & 2.2849(0.006)   \\
$\Sigma_c$ & $\{uu\}c$ & $1(\frac12^+)$ & 2.554 (0.45) & 2.439 & 2.452  & 2.4513(0.007)  \\
$\Xi_c$ & $[us]c$ & $\frac12(\frac12^+)$& 2.594 (0.45) &  2.481 & 2.473  & 2.4663(0.0014)  \\
$\Xi'_c$ & $\{us\}c$ & $\frac12(\frac12^+)$& 2.702 (0.44) & 2.578 &  2.599 & 2.5741(0.0033)\\
$\Omega_c$ &  $\{ss\}c$  & $0(\frac12^+)$& 2.829 (0.43) &  2.698 & 2.678 & 2.6975(0.0026)\\
$\Lambda_b$ & $[ud]b$ & $0(\frac12^+)$& 5.687 (0.45) &  5.622 & 5.672  & 5.6240(0.009)   \\
$\Sigma_b$ & $\{uu\}b$ & $1(\frac12^+)$ & 5.843 (0.42) & 5.805 & 5.847 &  \\
$\Xi_b$ & $[us]b$ & $\frac12(\frac12^+)$& 5.882 (0.41) &  5.812 &  5.788 &\\
$\Xi'_b$ & $\{us\}b$ & $\frac12(\frac12^+)$& 5.989 (0.40) & 5.937 &  5.936 &\\
$\Omega_b$ & $\{ss\}b$ & $0(\frac12^+)$& 6.115 (0.39) & 6.065 &  6.040 &\\
\bottomrule
\end{tabular}
\label{Table.baryon}
\end{table}

\subsection{Heavy tetraquarks}

In our calculations for heavy tetraquarks we have used the masses of diquark (antidiquark) and form factor parameters of Ref. \cite{ebert2006masses} which are given in Table \ref{Table.pot_parameters}.

Our numerical results for the masses of charm ($cq\bar{c}\bar{q}$ and $cs\bar{c}\bar{s}$) and bottom ($bq\bar b\bar{q}$ and $bs\bar b\bar{s}$) tetraquarks
for $s-$, $p-$ and $d-$wave channels with total spin ${\cal S}=0$ are listed in Tables \ref{Table.cqcq_cscs} and \ref{Table.bqbq_bsbs}. The tetraquark masses are calculated for scalar $S\bar{S}$ and axial-vector $A\bar{A}$ diquark--antidiquark contents.
We have also calculated the expectation value of the relative distance between $D\bar D$ pair which can provide an estimate of the size of the tetraquarks.

\begin{table}[hbt]
\caption{Masses of charm diquark--antidiquark states in units of GeV, calculated by nonrelativistic R--matrix method compared to nonrelativistic Lippmann--Schwinger (LS) calculations and relativistic EFG results. The numbers in parentheses are the expectation value of the relative distance between $D\bar D$ pair $\evdel{r}$ in units of fm.}
\centering
\begin{tabular}{ccccccccccccccc}
\toprule
state & \multicolumn{4}{c}{$cq\bar{c}\bar{q}\,(S\bar{S})$} & & \multicolumn{4}{c}{$cq\bar{c}\bar{q}\,(A\bar{A})$} \\ \cline{2-5} \cline{7-10}
 & R--matrix & LS \cite{Hadizadeh-PLB753} & EFG~\cite{ebert2006masses,Ebert_EPJC58} & LS \cite{Hadizadeh-EPJC75} && R--matrix & LS \cite{Hadizadeh-PLB753} & EFG~\cite{ebert2006masses,Ebert_EPJC58}  & LS \cite{Hadizadeh-EPJC75} \\ \hline
$1s$ & 3.885 (0.35) & 3.792 & 3.812 & 3.885 && 4.013 (0.35) & 3.919 & 3.852 & 4.013 \\
$1p$  & 4.268 (0.55) & 4.262 & 4.244 & && 4.388 (0.54) & 4.374 &4.350 \\
$2s$  & 4.461 (0.70) & 4.419 & 4.375 & && 4.580 (0.69) & 4.535 & 4.434 \\
$1d$  & 4.553 (0.73) & 4.556 & 4.506 & && 4.669 (0.72) & 4.668 & 4.617 \\
$2p$  & 4.708 (0.85) & 4.697 & 4.666 & && 4.823 (0.84) & 4.816 & 4.765 \\
$3s$  & 4.873 (0.99) & 4.843  & & && 4.988 (0.97) & 4.944  \\
$2d$  & 4.932 (1.00) & 4.933 & & && 5.044 (0.99) & 5.037  \\
$3p$  & 5.073 (1.10) & 5.062 & & && 5.184 (1.09) & 5.184  \\
\hline
state & \multicolumn{4}{c}{$cs\bar{c}\bar{s}\,(S\bar{S})$} & & \multicolumn{4}{c}{$cs\bar{c}\bar{s}\,(A\bar{A})$} \\ \cline{2-5} \cline{7-10}
 & R--matrix & LS \cite{Hadizadeh-PLB753} & EFG~\cite{ebert2006masses,Ebert_EPJC58} & LS \cite{Hadizadeh-EPJC75}  && R--matrix & LS \cite{Hadizadeh-PLB753} & EFG \cite{ebert2006masses,Ebert_EPJC58} & LS \cite{Hadizadeh-EPJC75}  \\ \hline
 $1s$  & 4.117 (0.35) & 4.011  & 4.051 & 4.117 && 4.249 (0.34) & 4.139  & 4.110 & 4.250 \\
$1p$  & 4.490 (0.54) & 4.490  & 4.466 & && 4.617 (0.53) & 4.616 & 4.582  \\
$2s$  & 4.681 (0.68) & 4.620  & 4.604 & && 4.808 (0.68)  & 4.744 & 4.680 \\
$1d$  & 4.770 (0.71) & 4.770 & 4.728 & && 4.893 (0.71) & 4.894 & 4.847 \\
$2p$  & 4.922 (0.83) & 4.920  & 4.884 & && 5.045 (0.82) & 5.041  & 4.991 \\
$3s$  & 5.086 (0.95) & 5.039 & &&& 5.208 (0.95) & 5.160 \\
$2d$  & 5.142 (0.98) & 5.143 & &&& 5.262 (0.97) & 5.263  \\
$3p$  & 5.281 (1.08) & 5.276 & &&& 5.399 (1.07) & 5.394  \\
\bottomrule
\end{tabular}
\label{Table.cqcq_cscs}
\end{table}

\begin{table}[hbt]
\caption{The same as Table \ref{Table.cqcq_cscs}, but for bottom tetraquarks.}
\centering
\begin{tabular}{cccccccccccccccc}
\toprule
state & \multicolumn{3}{c}{$bq\bar b \bar q\,(S\bar{S})$} & & \multicolumn{3}{c}{$bq\bar b \bar q\,(A\bar{A})$} \\ \cline{2-4} \cline{6-8}
 & R--matrix & LS \cite{Hadizadeh-PLB753} & EFG~\cite{ebert2006masses,Ebert_MPLA24} & & R--matrix & LS \cite{Hadizadeh-PLB753} & EFG~\cite{ebert2006masses,Ebert_MPLA24} \\ \hline
$1s$  & 10.482 (0.23) & 10.426 & 10.471 & & 10.527 (0.23) & 10.469  & 10.473 \\
$1p$  & 10.814 (0.38) & 10.813  & 10.807  && 10.858 (0.38) & 10.856  & 10.850 \\
$2s$  & 10.942 (0.47) & 10.914  & 10.917  && 10.986 (0.48) & 10.958  & 10.942 \\
$1d$  & 11.034 (0.51) & 11.034 &  11.021 && 11.077 (0.51) & 11.077 &  11.064  \\
$2p$  & 11.142 (0.59) & 11.140  & 11.122 && 11.185 (0.59) & 11.183 & 11.163 \\
$3s$  & 11.252 (0.68) & 11.230 & && 11.295 (0.68) & 11.273  \\
$2d$  & 11.311 (0.70) & 11.310 & && 11.354 (0.70) & 11.354  \\
$3p$  & 11.409 (0.78) & 11.406 & && 11.452 (0.78) & 11.450   \\
\hline
state & \multicolumn{3}{c}{$bs\bar b \bar s\,(S\bar{S})$} & & \multicolumn{3}{c}{$bs\bar b \bar s\,(A\bar{A})$} \\ \cline{2-4} \cline{6-8}
 & R--matrix & LS \cite{Hadizadeh-PLB753} & EFG~\cite{ebert2006masses,Ebert_MPLA24} & & R--matrix & LS \cite{Hadizadeh-PLB753} & EFG~\cite{ebert2006masses,Ebert_MPLA24} \\ \hline
 $1s$ & 10.691 (0.23) & 10.629 & 10.662  && 10.732 (0.23) & 10.668  & 10.671 \\
$1p$ & 11.017 (0.38) & 11.015 & 11.002  && 11.056 (0.38) & 11.054  & 11.039  \\
$2s$ & 11.146 (0.48) & 11.116 & 11.111 &&  11.186 (0.48) & 11.155  & 11.133  \\
$1d$ & 11.235 (0.51) & 11.235 &  11.216  && 11.275 (0.51) & 11.274 & 11.255  \\
$2p$ & 11.343 (0.59) & 11.340  & 11.316  && 11.382 (0.59) & 11.379 & 11.353  \\
$3s$ & 11.454 (0.68) & 11.430 & && 11.493 (0.68) & 11.469  \\
$2d$ & 11.511 (0.70) & 11.511 & && 11.550 (0.70) & 11.549  \\
$3p$  & 11.608 (0.78) & 11.606 & && 11.647 (0.77) & 11.645 \\
\bottomrule
\end{tabular}
\label{Table.bqbq_bsbs}
\end{table}

We have compared our results for tetraquark masses with recent results obtained in momentum and configuration spaces by solution of the nonrelativistic Lippmann--Schwinger integral equation \cite{Hadizadeh-PLB753,Hadizadeh-EPJC75} and also with those of previous relativistic studies by EFG reported in Refs. \cite{ebert2006masses,Ebert_EPJC58,Ebert_MPLA24}.

It indicates that the calculated masses for the ground state of charm tetraquarks (i.e. $cq\bar{c}\bar{q}$ and $cs\bar{c}\bar{s}$) with corresponding results from LS (in configuration space) \cite{Hadizadeh-EPJC75} shows that they are in excellent agreement.
Our results are also in good agreement with those of LS (in momentum space) and EFG with a relative percentage difference estimated to be at most $2.5 \,\%$ and $4\,\%$, respectively.
While our R--matrix calculations and LS results are both done in a nonrelativistic spin--independent scheme, there is some difference between the results obtained in momentum and configuration spaces.
Clearly, the difference between our R--matrix results in configuration space and EFG results in momentum space is larger, and comes from the relativistic effects and also spin contribution in the $D\bar{D}$ interaction which appears in spin--orbit, spin--spin and tensor spin--space terms \cite{Ebert_EPJC58}.
As we have shown in Ref. \cite{Hadizadeh-PLB753} the relativistic effect leads to a small reduction in the mass of heavy tetraquarks and decreases the masses of charm and bottom tetraquarks by less than $2\,\%$ and $0.2\,\%$, respectively, whereas the spin contribution may lead to small decrese or increase in the masses of tetraquarks.

\begin{figure}[hbt]
 \caption{The $x-z$ cross section of the $s-$, $p-$ and $d-$wave probability densities $|\Psi_{nlm}({\mathbf r})|^2$ for $bq\bar b\bar{q}$ tetraquark in $S\bar{S}$ state for $m=0$.}
 \centering
 \includegraphics[width=.8\textwidth]{bqbq_SS_density_group.pdf}
 \label{Fig.densities}
\end{figure}

In Fig. \ref{Fig.densities} we have shown few examples of the probability density $|\Psi_{nlm}({\mathbf r})|^2$ of $bq\bar b\bar{q}$ tetraquark in $S\bar{S}$ state for $s$, $p$ and $d$ channels.
As we can see the tetraquark probability functions for higher states have been expanded to larger distances which leads to larger expectation value for the relative distance between $D \bar D$ pair.
As we can see in Tables \ref{Table.cqcq_cscs} and \ref{Table.bqbq_bsbs}, the expectation value of the relative distance between $D \bar D$ pair is almost the same for scalar $S\bar{S}$ and axial-vector $A\bar{A}$ diquark--antidiquark contents.
It is larger for higher states and its size changes roughly with a factor of 3 from $1s$ to $3p$ state.
In Table \ref{Table.cqcq_experiment}, we have compared our results for the masses of charm and bottom tetraquarks with the possible experimental candidates. They are in good agreement with a relative difference below $3.4\,\%$.

\begin{table}[hbt]
\caption{Comparison of our numerical results for the masses of charm and bottom diquark--antidiquark states, calculated by R--matrix method and possible experimental candidates. The expectation value $\evdel{r}$ for diquark--antidiquark pair distance is also calculated. }
\centering
\begin{tabular}{cccccccccccccc}
\toprule
& \multirow{2}{*}{state} & \multicolumn{2}{c}{R--matrix Theory} & & \multicolumn{2}{c}{Experiment} \\ \cline{3-4} \cline{6-7}
&  & Mass (MeV) & $\evdel{r}\, (\text{fm})$ & & Exp. candidate & Mass (MeV)     \\ \hline
 \multirow{9}{*}{\rotatebox[origin=c]{90}{$cq\bar{c}\bar{q}\,(S\bar{S})$}} \\
& & & & & &
 \multirow{3}{*}{
$
\left\{
\begin{array}{l l}
 4259 \pm 8^{+2}_{-6} & \text{\cite{PhysRevLett.95.142001}} \\
 4247 \pm 12^{+17}_{-32} & \text{\cite{PhysRevLett.99.182004}} \\
 4284 ^{+17}_{-16} \pm 4 & \text{\cite{He-PRD74}} \\
\end{array}
\right.
$
} \\
& $1p$ & 4268 & 0.5500 & & $Y(4260)$ &\\
& & & & & \\
\\
& & & & & &
 \multirow{3}{*}{
 $
\left\{
\begin{array}{l l}
 4664\pm11\pm5 & \text{\cite{PhysRevLett.99.142002}} \\\\
 4634^{+8+5}_{-7-8} & \text{\cite{PhysRevLett.101.172001}}
\end{array}
\right.
$
} \\
& $2p$ & 4708 &0.8491 & & $Y(4660)$ &\\
& & & & & \\ \\
\hline
 \multirow{9}{*}{\rotatebox[origin=c]{90}{$cq\bar{c}\bar{q}\,(A\bar{A})$}} \\
& & & & & &
 \multirow{3}{*}{$
\left\{
\begin{array}{l l}
 4361\pm9\pm9 & \text{\cite{PhysRevLett.99.142002}} \\
 4324\pm24 & \text{\cite{PhysRevLett.98.212001}} \\
 4355^{+9}_{-10} \pm 9 & \text{\cite{{Liu_PRD78}}}
\end{array}
\right.
$} \\
& $1p$ & 4388 & 0.5449 && $Y(4360)$ &\\
& & & & & \\
\\
& & & & & &
 \multirow{3}{*}{$
\begin{array}{l l}
 4433\pm4\pm2 & \text{\cite{PhysRevLett.100.142001} }
\end{array}
$} \\
& $2s$ & 4580 &0.6924 && $Z(4430)$ &\\
& & & & & \\
\hline
 \multirow{9}{*}{\rotatebox[origin=c]{90}{$bq\bar b\bar{q}\,(S\bar{S})$}} \\
 & & & & & &
 \multirow{3}{*}{$
\left\{
\begin{array}{l l}
 10876 \pm 2  & \text{\cite{Aubert_PRL102}} \\ \\
 10865 \pm 8 & \text{\cite{Yao_JPG33}}
\end{array}
\right.
$}
\\
& $1p$ & 10814 & 0.3780 && $Y(10860)$   \\  \\ \\
 & & & & & &
 \multirow{3}{*}{$
\left\{
\begin{array}{l l}
 10996 \pm 2  & \text{\cite{Aubert_PRL102}} \\ \\
 11019 \pm 8 & \text{\cite{Yao_JPG33}}
\end{array}
\right.
$}
\\
& $2p$ & 11142 & 0.5948 && $Y(11020)$   \\ \\
\bottomrule
\end{tabular}
\label{Table.cqcq_experiment}
\end{table}

\subsection{Pentaquarks}

By successful application of the R--matrix method for diquark, triquarks and tetraquarks we have also implemented it to study the pentaquarks as bound states of diquark--antitriquark systems (see Fig. \ref{Fig.pentaquark}).

In this study we have used two models of $D\bar T$ interaction. In the following sections we present the models and our numerical results.

\subsubsection{One--pion exchange potential}

One--pion exchange potential (OPEP) acting between a nucleon and a heavy meson ($D$ or $B$) given as \cite{Cohen_PRD72}
\begin{eqnarray}
V_{\pi} (r) =\left\{
                \begin{array}{ll}
                 V_0 ,  & r< r_0
                  \\ \\
                 2 {\bf I}_N \cdot {\bf I}_H  \biggl ( S_{12} V_T(r) + 2 {\bf S}_N \cdot {\bf S}_l V_c(r)  \biggr ) , & r > r_0
                \end{array}
              \right.
\end{eqnarray}
where
\begin{eqnarray}
V_0 &=&  -62.79 \quad \text{or} \quad  -276 \, \text{MeV} , \cr
V_T(r) &=& \dfrac{g_A g_H m_{\pi}^2 }{2 \pi f_{\pi}^2}  \dfrac{e^{-m_{\pi}r}}{6r}  \left ( \dfrac{3}{m_{\pi}^2 r^2} + \dfrac{3}{m_{\pi}r} +1 \right ), \cr
V_c(r) &=& \dfrac{g_A g_H m_{\pi}^2 }{2 \pi f_{\pi}^2}  \dfrac{e^{-m_{\pi}r}}{3r}  ,
\end{eqnarray}
and the labels are:
nucleon isospin $I_N$, heavy-meson isospin $I_H$ (total isospin ${\bf I} = {\bf I}_N + {\bf I}_H$), tensor force $S_{12}$, tensor potential $V_T (r)$, nucleon spin $S_N$, light quark in heavy meson spin $S_l$ (sum of nucleon spin and light quark spin ${\bf K} = {\bf S}_N + {\bf S}_l$) and central potential $V_c(r)$.
The parameters of OPEP potential are given in Table \ref{Table.OPEP}.

\begin{figure}[hbt]
 \caption{Pentaquark system as a diquark--antitriquark bound state.}
\centering
 \includegraphics[width=.65\textwidth]{pentaquark}
 \label{Fig.pentaquark}
\end{figure}

\begin{table}[hbt]
\caption{The parameters of OPEP potential used in the heavy pentaquark calculations.}
\centering
\begin{tabular}{lllllllll}
\toprule
parameter &&& value \cr
\hline
$g_A$ &&& $1.27$ \cr
$f_{\pi}$ &&& $131$ MeV \cr
$g_H$ &&& $+0.59 $ \cr
$m_{\pi}$ &&& $138$ MeV \cr
$m_N$ &&& $938.92$ MeV \cr
$m_B$ &&& $5279$ MeV \cr
$m_D$ &&& $1867$ MeV \cr
\bottomrule
\end{tabular}
\label{Table.OPEP}
\end{table}

Our results for the masses of the pentaquarks composed of a nucleon and a meson ($B$ and $D$ mesons) for $l=0$, $S=\frac{3}{2}$ and $J^P=\frac{3}{2}^+$ channel are given in Table \ref{tab:B-D_meson}. As we can see our results with R--matrix method are in excellent agreement with Cohen et al. results \cite{Cohen_PRD72}.

\begin{table}[hbt]
\caption{The masses of pentaquarks composed of a nucleon and a meson for $l=0$, $S=\frac{3}{2}$ and $J^P=\frac{3}{2}^+$ channel. Upper panel: nucleon + $B$ meson, lower panel: nucleon plus $D$ meson. All masses are in GeV\@.  Column
A: constant potential,
$V_0=-276$~MeV and $r_0 = 1$~fm;
B: constant potential,
$V_0=-62.79$~MeV and $r_0 = 1.5$~fm.}
\centering
\begin{tabular}{lllllllll}
\toprule
& \multicolumn{2}{c}{A} && \multicolumn{2}{c}{B} \\
\cmidrule{2-3} \cmidrule{5-6}
 & \cite{Cohen_PRD72} & R--matrix && \cite{Cohen_PRD72} & R--matrix  \\
 \hline
 $I=0$ &  6.077  & 6.076  && 6.202 & 6.201  \\
 $I=1$ &  6.077 & 6.076 && 6.203 & 6.202  \\
 \hline
 $I=0$ & 2.689 & 2.688 && 2.797 & 2.795   \\
$I=1$ & 2.691 & 2.689  && 2.797 & 2.796  \\
\bottomrule
\end{tabular}
\label{tab:B-D_meson}
\end{table}

\subsubsection{Cornell Potential}
For the second test of Pentaquark calculations, the nonrelativistic linear--plus--Coulomb Cornell potential is used. It has the following form \cite{Barnes_PRD72}

\begin{equation} \label{Cornell_penta}
V(r) = V_c(r) + V_{spin-spin}+ V_{spin-orbit} + V_{tensor}
\end{equation}
where
\begin{eqnarray}
V_c(r) &=& -\dfrac{4}{3} \dfrac{\alpha_s}{r} + br \cr
V_{spin-spin} &=& \dfrac{32 \pi \alpha_s}{9 m_D m_T} \left ( \dfrac{\sigma}{\sqrt{\pi}}  \right )^3  e^{-\sigma^2 r^2} {\bf S}_D \cdot  {\bf S}_{\bar T} \cr
V_{spin-orbit} &=& \dfrac{1}{m_D m_T} \left (  \dfrac{2\alpha_s}{r^3} - \dfrac{b}{2r} \right ) {\bf L} \cdot  {\bf S} \cr
V_{tensor} &=& \dfrac{1}{m_D m_T} \dfrac{4\alpha_s}{r^3} T,
\end{eqnarray}
and spin--spin, spin--orbit and tensor operators can be calculated as
\begin{eqnarray}
{\bf S}_D \cdot  {\bf S}_{\bar T} &=& \dfrac{1}{2} \left ( {\bf S}^2 - {\bf S}_D^2 - {\bf S}_{\bar T}^2  \right ), \cr
{\bf L} \cdot  {\bf S} &=& \dfrac{1}{2} \left ( {\bf J}^2 - {\bf L}^2 - {\bf S}^2  \right ), \cr
\langle ^3L_J | T | ^3L_J \rangle &=&
 \left\{
                \begin{array}{ll}
                  \dfrac{-L}{6(2L+3)},  & J=L+1
                  \\
                 \dfrac{1}{6}, & J=L
                  \\
                  -\dfrac{L+1}{6(2L-1)},  & J=L-1.
                \end{array}
              \right.
\end{eqnarray}
In Table \ref{Table.cornel-pentaquark}, the parameter of Cornell potential used in pentaquark calculations are given.

\begin{table}[hbt]
\caption{The parameters of Cornell potential used in heavy pentaquark calculations.}
\centering
\begin{tabular}{lllllllll}
\toprule
parameter &&& value \cr
\hline
$\alpha_s$ &&& $0.5461$ \cr
$b$ &&& $0.1425$ GeV$^2$ \cr
$\sigma$ &&& $1.0946$ GeV \cr
$m_D$ &&& $1.860$ GeV \cr
$m_{\bar T}$ &&& $2.286$ GeV \cr
\bottomrule
\end{tabular}
\label{Table.cornel-pentaquark}
\end{table}

Our results for the masses of charmoniumlike pentaquark $P_c^+$ in $s-$ and $p-$wave channels are given in Table \ref{Table.Pc+}.
In our calculations we have ignored tensor force and we have considered central, spin--spin and spin--orbit terms of diquark--antitriquark interaction of Eq. (\ref{Cornell_penta}).
Beside the small difference between our results for the masses of $s-$ and $p-$wave charmoniumlike pentaquarks and Lebed's results \cite{Lebed_PLB749}, which is about $5 \%$ and comes from neglected tensor force, the diquark--antitriquark separation $\evdel{r}$ for $s-$wave channel with value of $0.306$ fm is almost half of the obtained separation by Lebed with value of $0.64$ fm. For the $p-$wave channel, neglecting the tensor force leads to the relative difference of $14\%$ in calculated separations.

\begin{table*}[hbt]
\caption{Masses of charmoniumlike pentaquark $P_c^+$ in a diquark--antitriquark picture in units of MeV, calculated by nonrelativistic R--matrix method compared to Lebed's results. The diquark--antitriquark separation $\evdel{r}$ is also calculated in units of fm.}
\centering
\begin{tabular*}{\textwidth}{@{\extracolsep{\fill}}lllllllllllllllll@{}}
\toprule
state & $J^P$ & potential & Mass (MeV) & $\evdel{r}$ (fm)   \\ \hline
$s-$wave   & $\frac{3}{2}^-$  \\
& & & \multicolumn{2}{c}{R--matrix} \\  \cmidrule{4-5}
                  &  & $V_c$  & 4112 &  0.292 \\
                  &  & $V_c + V_{spin-spin} (+ V_{spin-orbit} ) $  & 4151 &  0.306  \\
 &  &  &  \multicolumn{2}{c}{Theory \cite{Lebed_PLB749}} \\ \cmidrule{4-5}
  &  & $V(r)$ & $4380$ & $0.64$  \\
   &  &  &  \multicolumn{2}{c}{\cite{Aaij_PRL115}} \\ \cmidrule{4-5}
  &  & Experimental candidate $P_c^+(4380)$ & $4380\pm8\pm29$  \\
  \hline
$p-$wave   & $\frac{5}{2}^-$   \\
& & & \multicolumn{2}{c}{R--matrix} \\  \cmidrule{4-5}
                  &  & $V_c$  & 4593 & 0.554  \\
                  &  & $V_c + V_{spin-spin}$  & 4597 & 0.559   \\
                  &  & $V_c + V_{spin-spin} + V_{spin-orbit} $ & 4633 &  0.600   \\
                   &  &  &  \multicolumn{2}{c}{Theory \cite{Lebed_PLB749}} \\ \cmidrule{4-5}
  &  & $V(r)$ & $4450$ & $0.70$  \\
     &  &  &  \multicolumn{2}{c}{\cite{Aaij_PRL115}} \\ \cmidrule{4-5}
  &  & Experimental candidate $P_c^+(4450)$ & $4449.8\pm1.7\pm2.5$  \\
\bottomrule
\end{tabular*}
\label{Table.Pc+}
\end{table*}

In conclusion, we have implemented the R--matrix method to calculate the mass spectra of heavy quarkonia, baryons, tetraquarks and pentaquarks in the two-body picture.
Our numerical results for the masses of heavy charm and bottom few--quarks even by neglecting the relativistic effects, are in good agreement with other theoretical predictions and also with available experimental data.

\begin{acknowledgements}
We would like to thank Y. Suzuki for helpful discussions about the application of R--matrix method to two--body problems, to R. F. Lebed for helpful comments on Cornel potential used in Pentaquark calculations and to Reza Abedi, 3D modeler and character artist, for creating the pentaquark picture.
M.A.S. acknowledges the financial support from the Brazilian agency CAPES and M.R.H. acknowledges the partial support by National Science Foundation under Contract No. NSF-HRD-1436702 with Central State University and by the Institute of Nuclear and Particle Physics at Ohio University.

\end{acknowledgements}