# Improving Quantum Computation
by Optimized Qubit Routing

Friedrich Wagner Fraunhofer Institute for Integrated Circuits IIS, Erlangen,
Germany<EMAIL_ADDRESS>Andreas Bärmann Department of Data
Science, University of Erlangen-Nuremberg, Germany<EMAIL_ADDRESS>Frauke Liers Department of Data Science, University of Erlangen-Nuremberg,
Germany<EMAIL_ADDRESS>Markus Weissenbäck Fraunhofer Institute for
Integrated Circuits IIS, Erlangen, Germany,
<EMAIL_ADDRESS>
###### Abstract

In this work we propose a high-quality decomposition approach for qubit
routing by swap insertion. This optimization problem arises in the context of
compiling quantum algorithms formulated in the circuit model of computation
onto specific quantum hardware. Our approach decomposes the routing problem
into an allocation subproblem and a set of token swapping problems. This
allows us to tackle the allocation part and the token swapping part
separately. Extracting the allocation part from the qubit routing model of
Nannicini et al. ([19]), we formulate the allocation subproblem as a binary
linear program. Herein, we employ a cost function that is a lower bound on the
overall routing problem objective. We strengthen the linear relaxation by
novel valid inequalities. For the token swapping part we develop an exact
branch-and-bound algorithm. In this context, we improve upon known lower
bounds on the token swapping problem. Furthermore, we enhance an existing
approximation algorithm which runs much faster than the exact approach and
typically is able to determine solutions close to the optimum. We present
numerical results for the fully integrated allocation and token swapping
problem. Obtained solutions may not be globally optimal due to the
decomposition and the usage of an approximation algorithm. However, the
solutions are obtained fast and are typically close to optimal. In addition,
there is a significant reduction in the number of artificial gates and output
circuit depth when compared to various state-of-the-art heuristics. Reducing
these figures is crucial for minimizing noise when running quantum algorithms
on near-term hardware. As a consequence, using the novel decomposition
approach leads to compiled algorithms with improved quality. Indeed, when
compiled with the novel routing procedure and executed on real hardware, our
experimental results for quantum approximate optimization algorithms show an
significant increase in solution quality in comparison to standard routing
methods.
Keywords: Combinatorial Optimization, Integer Programming, Approximation
Algorithms, Quantum Compilation

##  1 Introduction

_Qubit routing by swap insertion_ is a subroutine in the process of compiling
quantum algorithms onto specific hardware. Quantum algorithms are usually
formulated in the circuit model of quantum computation, see e.g. [20] for an
introduction. Such a circuit consists of wires representing qubits as well as
gates representing operations applied to subsets of qubits. We refer to the
qubits in the circuit as _logical_ qubits, in contrast to the _physical_
qubits in quantum hardware. Furthermore, we assume that the circuit only
contains gates acting on at most two qubits, i.e. single and two-qubit gates
(also known as “$\mathrm{SU}(4)$ circuits”). This might already be the case
for the input circuit to be compiled; otherwise this can be achieved using the
Solovay-Kitaev Theorem, see e.g. [20, pp. 188–202]. In the circuit model,
gates may be applied on any pair of logical qubits. However, this not the case
in currently available quantum processors. Here, two-qubit gates can only be
applied on specific pairs of physical qubits defined by the _hardware
connectivity graph_. Now the task is to choose an initial allocation of
logical qubits in the quantum circuit to physical qubits in the hardware graph
such that for each two-qubit gate the corresponding logical qubits are located
at neighbouring physical qubits. The example in Fig. 1 shows that this is not
always possible directly.

(a)

(b)

Figure 1: Example for a quantum circuit (a) consisting of four logical qubits
as well as single and two-qubit gates. In the context of this work, the
meaning of the gates is irrelevant. Only the subset of qubits involved in a
gate matters. The circuit is not compatible with the hardware connectivity
graph shown in (b).

In such a case, it is necessary to insert additional _swap gates_ into the
algorithm. They effectively change the allocation of logical qubits to
physical qubits. The result is an $\mathrm{SU}(4)$ circuit meeting the
connectivity restrictions which is equivalent to the original one when taking
into account the permutation resulting from initial allocation and swap gates.
A simple example for this procedure is shown in Figure 2.

Figure 2: Example for routing by swap insertion. The input circuit and
hardware graph are taken from Figure 1. The input circuit in Figure 1a does
not meet the hardware connectivity requirement. The given output circuit is
equivalent but fulfils the hardware connectivity. The equivalence holds under
the initial allocation of logical qubits $\\{q_{0},q_{1},q_{2},q_{3}\\}$ to
physical qubits $\\{0,1,2,3\\}$ given by $q_{0}\mapsto 3,q_{1}\mapsto
2,q_{2}\mapsto 1,q_{3}\mapsto 0$. After the first swap the allocation changes
to $q_{0}\mapsto 3,q_{1}\mapsto 2,q_{2}\mapsto 0,q_{3}\mapsto 1$, and after
the second swap to $q_{0}\mapsto 3,q_{1}\mapsto 2,q_{2}\mapsto 1,q_{3}\mapsto
0$.

When solving the routing by swap insertion problem, two possible objective
functions arise naturally: either one aims to minimize the execution time on
the quantum processor, which can be achieved by minimizing the output circuit
depth, or one can minimize the total gate number to be executed, which amounts
to minimizing the number of swap gates added. Intuitively, both objectives are
correlated. However, it is easy to construct examples where optimal-depth
solutions are suboptimal in terms of gate count and vice versa. For an
experimental study of the influence of depth and gate count on the performance
of quantum algorithms, we refer to [19] and [12].

##### Existing methods for qubit routing.

Many compilation procedures based on routing by swap insertion have been
developed so far. Some of them may be classified by their local working
principle: they define an initial mapping of logical to physical qubits,
iterate through the circuit in temporal order and change the allocation by
adding swaps when gates are encountered which cannot be applied in the current
allocation. Of course, this general principle allows for many different
sophisticated procedures, in particular when it comes to choosing the initial
allocation. Cowtan et al.’s _Tket compiler_ (see [6, 23]) belongs to this
class. It uses a heuristic cost function to choose the next allocation.
Zuhlener et al. (see [30]) employ an A*-search algorithm to define the next
allocation. Li et al.’s _Sabre compiler_ (see [15]) uses the reversibility of
quantum circuits to perform a bidirectional search for a good initial mapping.
The default compiler in IBMs SDK _Qiskit_ , _Stochastic swap_ , based on work
by Bravi (see [2]), as well as the default compiler in Googles SDK _Circ_ ,
_Greedy router_ (see [8]) also belong to this class.

A connection from routing by swap insertion to the problems of token swapping
and subgraph isomorphism has already been established and employed by Siraichi
et al. in [22, 21] and by Childs et al. in [3]. In their work, the subgraph
isomorphism problem arises in the context of searching for high-quality qubit-
mappings, whereas routing between mappings is translated to a token swapping
problem.

Two approaches to routing by swap insertion that do not rely on local swap
insertion but consider the whole circuit instead are the SAT-based approach by
Wille et al. (see [27]) and the approach by Nannicini et al. based on an
integer programming model (see [19]). These approaches suffer from exponential
growth of problem size and NP-hard problem complexity. In practice, the
running times exceed reasonable limits already for moderate circuit sizes.

Furthermore, there are quantum compilation methods in the literature which do
not rely on swap insertion: Kissinger and Meijer-Van De Griend consider
circuits only built from CNOT gates and solve the compiling problem by finding
an equivalent hardware compatible CNOT circuit (see [14]). Moro et al. use
reinforcement learning for quantum compilation (see [18]).

Kim (see [13]) compares some of the aforementioned routing heuristics on
benchmarks with known optimal solution. His results reveal a significant
margin for improvement.

##### Our contribution.

We provide a high-quality solution method that efficiently exploits the
capabilities of integer programming to solve the routing by swap insertion
problem. The key is to decompose the problem into two separate subproblems of
which the first one is solved via integer programming. It determines an
allocation sequence that is compatible with the hardware constraints. Herein,
we employ a cost function that is a lower bound on the total number of swaps
required for routing. This subproblem is called _token allocation problem_.
For its solution, we develop a simplification of the binary model introduced
by Nannicini et al. in [19], which results in a significant reduction in
problem size. Once an optimal allocation sequence has been found, each pair of
subsequent allocations defines a token swapping instance. These token swapping
problems are solved by an approximation algorithm based on the work of Miltzow
et al. in [16]. Also, we develop an exact algorithm for token swapping. By
comparison of both solution techniques, we conclude that the approximation
algorithm typically delivers high-quality solutions in drastically reduced
solution time.

Re-targeting the gates according to the allocation sequence and inserting the
swap gates from the token swapping solutions finally yields the compiled
circuit. In [21] a similar decomposition approach is proposed which, however,
is based on dynamic programming. In contrast to other heuristic methods, we
are able to give bounds on the quality of the obtained solution. The numerical
results show that the proposed approach finds a close-to-optimal solution on
instances with known optimum and outperforms well-established heuristics. In
comparison to the exact approach from [19], on which the present work is
based, our method takes much less time. This makes the proposed routing
procedure applicable for practically relevant circuit sizes. Furthermore, our
experiments on actual quantum hardware show that quantum computation benefits
from the proposed routing method.

##### Structure.

This work is structured as follows. In Section 2 we introduce the token
allocation problem and the token swapping problem. The allocation problem is
further studied in Section 3. We show NP-hardness, introduce the binary model
and derive valid inequalities for its linear relaxation. In Section 4 we turn
to the token swapping problem. The efficient approximation algorithm employed
in our routing procedure as well as the exact branch and bound algorithm used
for benchmarking the former are developed. Section 5 discusses numerical
results for the token allocation model and the token swapping algorithm as
well as the integrated routing problem Additionally, we present experimental
results for example quantum algorithms executed on actual hardware when routed
with different methods. Finally, in Section 6 we give a conclusion and
indicate directions of further research.

##  2 Preliminaries and Definitions

We start by formally defining the token allocation problem as well as the
token swapping problem. While the latter is well known, the former has not
been extensively studied yet, to the authors’ best knowledge. It already
occurred in [19] and [22, 21] as part of a larger problem. In the context of
qubit routing, the token allocation problem can be informally described as
follows. Its input is a hardware graph representing the connectivity of the
physical qubits together with a quantum algorithm acting on a set of logical
qubits (“tokens”). The algorithm is described as a sequence of layers. In each
layer, a set of two-qubit gates needs to be performed in parallel. Now the
task is to find an allocation for each layer such that the gates can be
executed, i.e. logical qubits involved in gates are located at neighbouring
physical qubits.

We remark that grouping gates into layers needs to be performed by the routing
procedure, since an algorithm formulated as a circuit is not by itself divided
into layers but simply modeled as a sequence of individual gates. Grouping
gates into layers containing more than a single gate restricts the set of
feasible solutions, since changing allocation between gates in the same layer
is not allowed. However, grouping significantly reduces the number of layers
and thus the size of the token allocation problem.

For a given allocation sequence, we define its cost as the sum of all
distances logical qubits move on the hardware graph between subsequent
allocations. Herein, the distance a logical qubit moves between two
allocations is defined as the length of a shortest path connecting the
vertices in the hardware graph to which the qubit is allocated. The motivation
for this choice of objective is that the number of swaps needed for routing
between subsequent allocations $a,a^{\prime}\colon Q\xrightarrow{1:1}V$ is
bounded from below by

$\frac{1}{2}\sum_{q\in Q}d_{H}\left(a(q),a^{\prime}(q)\right),$ (1)

where $d_{H}\colon V\times V\rightarrow\mathbbm{Z}_{0}^{+}$ denotes the _edge
distance_ in $H$, i.e. the number of edges in a path of minimal length
connecting two nodes. A proof of this result on the token swapping problem is
found e.g. in [16]. This means that instead of minimizing the total number of
swaps required, we minimize a lower bound on this value which can be computed
much more easily.

We now define the token allocation problem formally as a combinatorial
optimization problem.

######  Definition 2.1 (Token allocation problem).

An instance of the _token allocation problem (TAP)_ is given by a tripel
$(H,Q,\Gamma)$, where

  * •

$H=(V,E)$ is a connected undirected graph,

  * •

$Q$ is a set of tokens of size $\lvert Q\rvert=\lvert V\rvert$, and

  * •

$\Gamma=G^{1},G^{2},\ldots,G^{L}$ is a finite sequence of sets of disjoint
token pairs, i.e.

    1. (i)

$G^{t}=\\{(p^{1}_{t},q^{1}_{t}),(p^{2}_{t},q^{2}_{t}),\ldots,(p^{\lvert
G^{t}\rvert}_{t},q^{\lvert G^{t}\rvert}_{t})\\}\subset Q\times
Q,\,q^{j}_{t}\neq p^{j}_{t}\quad\forall j\in\\{0,\ldots,\lvert G^{t}\rvert\\}$
and

    2. (ii)

$\\{p^{i}_{t},q^{i}_{t}\\}\cap\\{p^{j}_{t},q^{j}_{t}\\}=\emptyset\quad\forall
i\neq j\in\\{0,\ldots,\lvert G^{t}\rvert\\},\,\forall t\in\\{1,\ldots,L\\}$.

A feasible solution is given by a sequence $a_{1},\ldots,a_{L}$ of bijective
mappings $a_{t}\colon Q\xrightarrow{1:1}V$ such that at each time step $t$ the
token pairs in $G^{t}$ are allocated to neighbouring vertices in $H$:

$\displaystyle\forall 1\leq t\leq L:\,\forall(q_{1},q_{2})\in
G^{t}:\\{a_{t}(q_{1}),a_{t}(q_{2})\\}\in E.$

The cost of a feasible solution $a_{1},\ldots,a_{L}$ is defined as

$\frac{1}{2}\sum_{t=0}^{L-1}\sum_{q\in
Q}d_{H}\left(a_{t}(q),a_{t+1}(q)\right)\,.$ (2)

The goal is to find a feasible allocation minimizing this cost function.

In the context of qubit routing, the vertex set $V$ represents the physical
qubits while the edge set $E$ represents the set of physical qubit pairs
between which two-qubit gates can be applied. The logical qubits are
represented by the set $Q$, and the sets $G^{t}\in\Gamma$, $1\leq t\leq L$,
are the gates to be executed in parallel at time step $t$.

Next, we give an informal description of the token swapping problem. Its input
is an undirected, connected graph and a set of tokens. Initially, each vertex
holds a unique token. Further, each vertex is assigned a target token. Now the
task is to move tokens along the edges of the graph such that each vertex
holds its desired token. However, movement of tokens is restricted to swapping
adjacent tokens. Formally, the token swapping problem is defined as follows.

######  Definition 2.2 (Token swapping problem).

An instance of the _token swapping problem_ is given by a tupel
$(H,Q,a,a^{\prime})$, where

  * •

$H=(V,E)$ is a connected undirected graph,

  * •

$Q$ is a set of tokens with size $\lvert Q\rvert=\lvert V\rvert$,

  * •

$a\colon Q\xrightarrow{1:1}V$ is an initial allocation, and

  * •

$a^{\prime}\colon Q\xrightarrow{1:1}V$ is a final allocation.

For two nodes $i,j\in V$, $i\neq j$, the transposition $(i,j)$ is the unique
bijective mapping from $V$ onto itself interchanging exactly the two elements
$i$ and $j$. Let $\mathcal{T}_{H}$ denote the set of transpositions on $V$
restricted to $E$, i.e. $(i,j)\in\mathcal{T}_{H}\Leftrightarrow\\{i,j\\}\in
E$. A feasible solution $\mathcal{S}$ is a finite sequence of transpositions
(“swaps”) $\mathcal{S}=\tau_{1},\tau_{2}\ldots,\tau_{N}$, where
$\tau_{i}\in\mathcal{T}_{H}$ for $i\in\\{1,\ldots,N\\}$, such that

$a=\tau_{N}\circ\ldots\circ\tau_{2}\circ\tau_{1}\circ a^{\prime}\ .$ (3)

The cost of a feasible solution is given by its length $N$. The goal of the
token swapping problem is to find a feasible solution minimizing these costs.

There is a close link between Problems 2.1 and 2.2: every feasible solution to
a TAP instance defines a sequence of token swapping instances in a canonical
way. In the application to qubit routing, a solution to this sequence of
problems together with the TAP solution results in a solution to the qubit
routing by swap insertion problem. Note that this solution might not be
optimal in terms of used swaps.

##  3 An Analysis of the Token Allocation Problem

In the following, we analyse the TAP from Definition 2.1. After showing its
NP-hardness, we derive an integer programming formulation based on edge flows.
Additionally, we strengthen the model by introducing valid inequalities.

###  3.1 NP-hardness

To motivate the use of integer programming for the solution of the TAP, we
first show its NP-hardness. For this purpose, we construct a reduction of the
_subgraph isomorphism problem_ to the TAP. A similar reduction was already
sketched by Siraichi et al. in [22] to show NP-hardness of a problem called
_qubit assignment problem_. This problem is equivalent to the decision version
of the TAP. Here, we work out the reduction in detail.

######  Definition 3.1 (Subgraph Isomorphism Problem).

Given two graphs $H$, $H^{\prime}$, the _edge-induced (node-induced) subgraph
isomorphism problem (SGI)_ asks whether there is an edge-induced (node-
induced) subgraph of $H$ isomorphic to $H^{\prime}$.

The complexity of the node-induced SGI is well known.

######  Theorem 3.2 ([26]).

The node-induced subgraph isomorphism problem is NP-complete.

To transfer the NP-completeness to the edge-induced SGI, we need the notion of
a _line graph_.

######  Definition 3.3 (Line Graph).

The _line graph_ $\mathcal{L}(G)$ of a graph $G$ possesses a node for every
edge in $G$. Edges in $\mathcal{L}(G)$ connect two nodes if and only if the
corresponding edges in $G$ share a common node.

Due to space limitations we state the following result without proof.

######  Lemma 3.4.

Two graphs $H$ and $H^{\prime}$ are node-induced subgraph isomorphic if and
only if the corresponding line graphs $\mathcal{L}(H)$ and
$\mathcal{L}(H^{\prime})$ are edge-induced subgraph isomorphic.

Using this result, the node-induced SGI can be reduced to the edge-induced
SGI.

######  Corollary 3.5.

The edge-induced subgraph isomorphism problem is NP-complete.

For the reduction of edge induced SGI to TAP, we need

######  Definition 3.6 (Connectivity Graph).

For a finite set of tokens $Q$, let
$G\coloneqq\\{(q_{1},p_{1}),\ldots,(q_{n},p_{n})\\}\subset Q\times Q$, with
$q_{i}\neq p_{i}$ for $i=1,\ldots,n$ be a set of token pairs. Then the
_connectivity graph_ of $G$ is given by $C(G)=(\tilde{Q},E_{C})$ with

$\tilde{Q}\coloneqq\bigcup_{1\leq k\leq
n}\\{q_{k}\\}\cup\\{p_{k}\\}\quad\text{and}\quad
E_{C}\coloneqq\\{\\{q,p\\}\mid(q,p)\in G\\}.$ (4)

An example of a connectivity graph is depicted in Figure 3. The following
Lemma establishes a connection between the connectivity graph and the TAP. The
proof is not hard and omitted here.

Figure 3: Example for a connectivity graph. This example is constructed from
the quantum circuit in Figure 1. Logical qubits in the circuit are interpreted
as tokens, gates are interpreted as token pairs.

######  Lemma 3.7.

A TAP instance $(H,Q,\Gamma)$ has optimal value of $0$ if and only if there is
an edge-induced subgraph of $H$ isomorphic to the connectivity graph
$C(\Gamma)$.

With these preliminaries established, we now study the complexity of the TAP.

######  Theorem 3.8.

The token allocation problem is NP-hard.

###### Proof.

To show NP-hardness, we construct a polynomial reduction of the edge-induced
SGI to the decision version of the TAP. Let $H=(V,E)$ and
$H^{\prime}=(V^{\prime},E^{\prime})$ be two graphs. For the SGI instance
$(H,H^{\prime})$, we construct an TAP instance $(H,Q,\Gamma)$. Set $Q\coloneqq
V$ and construct $\Gamma$ as follows: choose an arbitrary order
$1,\ldots,\lvert E^{\prime}\rvert$ of the edges in $E^{\prime}$, and choose an
arbitrary order of the two vertices in each edge. This yields a sequence
$(p_{1},q_{1}),(p_{2},q_{2}),\ldots,(p_{\lvert E^{\prime}\rvert},q_{\lvert
E^{\prime}\rvert})$, $\\{p_{k},q_{k}\\}\in E^{\prime}$, $\forall
k\in\\{1,\ldots,\lvert E^{\prime}\rvert\\}$. Now, set
$G^{k}\coloneqq\\{(p_{k},q_{k})\\}$ for all $k\in\\{1,\ldots,\lvert
E^{\prime}\rvert\\}$ and $\Gamma\coloneqq G^{1},\ldots,G^{\lvert
E^{\prime}\rvert}$. Then we have $C(\Gamma)=H^{\prime}$. Finally, Lemma 3.7
implies that the TAP instance $(H,Q,\Gamma)$ has optimal value $0$ if and only
if there is an edge-induced subgraph isomorphism from $H^{\prime}$ to $H$. ∎

###  3.2 Network Flow Model

In [19], Nannicini et al. formulate routing by swap insertion as a network
flow problem. We adjust their model to the TAP. This results in a model
significantly reduced in size which, of course, is due to the fact that only a
subproblem of routing by swap insertion is modeled. Let $H=(V,E)$ be a
connected graph, and let $(H,Q,\Gamma)$ be a TAP instance. We introduce the
directed arc set $A_{H}\coloneqq\bigcup_{\\{u,v\\}\in
E}\\{(u,v)\\}\cup\\{(v,u)\\}$ as well as variables with the following
interpretations. The binary variables $x^{t}_{q,i,j}\in\\{0,1\\}$ take a value
of $1$ if qubit $q\in Q$ moves from node $i\in V$ to node $j\in V$ between
time steps $t$ and $t+1$ from $\\{1,\ldots,L\\}$, and $0$ otherwise. The
auxiliary binary variables $w_{q,i}^{t}\in\\{0,1\\}$ indicate by a value of
$1$ whether qubit $q\in Q$ is located at node $i\in V$ in time step
$t\in\\{1,\ldots,L\\}$, or $0$ if not. Further binary auxiliary variables
$y^{t}_{(p,q),(i,j)}\in\\{0,1\\}$ express whether gate $(p,q)\in G^{t}$ is
performed along edge $(i,j)$, value $1$, or not, value $0$. With these
notions, the TAP can be modeled by the following binary linear program: x∑t =
1L-1∑q∈Q∑i,j∈V×V ci,j xtq,i,j w^t_q,i =∑j∈Vxtq,i,j ∀1≤t < L,∀i ∈V,∀q ∈Q
w^t_q,i =∑j∈Vxt-1q,j,i ∀1 < t ≤L,∀i ∈V,∀q ∈Q ∑(i,j)∈AH yt(p,q),(i,j)=1
∀1≤t≤L,∀(p,q)∈G^t y^t_(p,q),(i,j)= w_p,i^t⋅w_q,j^t
∀1≤t≤L,∀(p,q)∈G^t,∀(i,j)∈A_H ∑i∈V wq,it= 1 ∀1≤t≤L,∀q∈Q ∑q∈Q wq,it= 1
∀1≤t≤L,∀i∈V w_q,i^t∈{0,1} ∀1≤t≤L,∀q∈Q,∀i∈V x^t_q,i,j∈{0,1}
∀1≤t≤L-1,∀q∈Q,∀(i,j)∈V×V y^t_(p,q),(i,j)∈{0,1} ∀1≤t≤L,∀(p,q)∈G^t,∀(i,j)∈A_H.
Constraints (3.2) and (3.2) ensure logical qubit conservation. Via Constraints
(3.2), we enforce that every gate is implemented. Constraints (3.2) demand
that a gate be implemented along an arc if and only if logical qubits are
located at the physical qubits of the arc. Finally, Constraints (3.2) and
(3.2) establish bijective mappings between logical and physical qubits while
Constraints (3.2),(3.2),(3.2) define the domains of the variables.

Model (3.2) can be illustrated by a time-expanded hardware graph $G$ which
contains a copy of the nodes in $H$ for every time step as well as edges
connecting subsequent time steps. An example is shown in Figure 4.

Figure 4: Example for the time-expanded hardware graph $G$ on which the
network flow formulation of the TAP is based. The instance is taken from the
example in Figure 1. Accordingly, there are $\lvert Q\rvert=4$ logical qubits
and $L=3$ time steps. The first number in a node labels the time step while
the second number labels the physical qubit.

For our further analysis of the model, we linearize Constraints (3.2) by
replacing them with

$\displaystyle y^{t}_{(p,q),(i,j)}$ $\displaystyle\leq w_{p,i}^{t}$
$\displaystyle y^{t}_{(p,q),(i,j)}$ $\displaystyle\leq w_{q,j}^{t}$ (5)
$\displaystyle y^{t}_{(p,q),(i,j)}$ $\displaystyle\geq
w_{p,i}^{t}+w_{q,j}^{t}-1.$ (6)

##### Optimal solution to the linear relaxation.

To illustrate the need for the introduction of strengthening inequalities, we
show that the linear programming relaxation of Model (3.2) is very weak in the
sense that for _any_ TAP instance there is an optimal LP solution with value
$0$. To see this, we assign the flow value $1/\lvert Q\rvert$ to all paths in
$G$ with zero costs. More formally, we set

$\displaystyle x^{t}_{q,i,j}$ $\displaystyle=\begin{cases}1/|Q|,&i=j\\\
0,&i\neq j\end{cases}\qquad\forall i,j\in V,\,\forall q\in Q,\,\forall 1\leq
t\leq L-1,$ $\displaystyle w^{t}_{q,i}$ $\displaystyle=1/\lvert
Q\rvert\quad\forall i\in V,\,\forall q\in Q,\,\forall 1\leq t\leq L$ and
$\displaystyle y^{t}_{(p,q),(i,j)}$ $\displaystyle=1/\lvert
A_{H}\rvert\quad\forall 1\leq t\leq L,\,\forall(p,q)\in
G^{t},\,\forall(i,j)\in A_{H}.$

It can be checked that this fully symmetrical, fractional solution is valid
for the linear relaxation of Model (3.2). Concerning Constraints (3.2) note
that $\lvert A_{H}\rvert\geq 2(\lvert Q\rvert-1)\geq\lvert Q\rvert$ and
$2/\lvert Q\rvert-1\leq 0$ hold, since $H$ is connected and we have $\lvert
Q\rvert\geq 2$.

###  3.3 Subgraph Isomorphism Inequalities

We will now derive valid inequalities for Model (3.2) to strengthen its LP
relaxation. To this end, we extend the result of Lemma 3.7, which we also used
in the proof of Theorem 3.8. For a given TAP instance $(H,Q,\Gamma)$, this
will enable us to identify a qubit subset $\tilde{Q}\subseteq Q$, an integer
distance $d\geq 0$ in $H$ as well as time steps $t_{0}$ and $t_{1}$ between
which at least one qubit pair in $\tilde{Q}$ has to move a total distance of
at least $d+1$. For the derivation we need

######  Definition 3.9 ($d$-th Relaxed Graph).

For an undirected, connected graph $H=(V,E)$ and an integer $d\geq 0$, we
define the _$d$ -th relaxed graph_ $H^{d}=(V,E^{d})$ as the graph defined by

$\displaystyle\\{i,j\\}\in E^{d}\iff d_{H}(i,j)\leq d+1\ .$

Note that $H^{0}=H$. An example for Definition 3.9 is shown in Figure 5.

(a) $H^{0}=H$

(b) $H^{1}$

(c) $H^{2}$

Figure 5: Example for the $d$-th relaxed graphs $H^{d}$ from Definition 3.9.
The graph $H$ is the four-qubit line graph from the example in Figure 1b.

We are now able to sketch the main idea for deriving valid inequalities.
Choose a starting time $t_{0}$ and an end time $t_{1}$ with $1\leq
t_{0}<t_{1}\leq L$ as well as a subset $\tilde{G}$ of the gates between
$t_{0}$ and $t_{1}$. Now, consider the connectivity graph
$C(\tilde{G})=(\tilde{Q},E_{C})$. We will relate it to the $d$-th relaxed
graph $H^{d}$ of $H$ for some integer $d\geq 0$. Namely, if there is no
isomorphism from $C(\tilde{G})$ to some edge-induced subgraph of $H^{d}$, we
know that there is no allocation such that all qubits involved in the gates of
$\tilde{G}$ are at most an edge distance of $d+1$ apart. Consequently, there
has to be at least one gate in $\tilde{G}$ such that the involved qubits are
more than $d+1$ apart at $t_{0}$. Between time steps $t_{0}$ and $t_{1}$, this
qubit pair will move a total distance of at least $d+1$. This leads us to the
following

######  Lemma 3.10 (Subgraph Isomorphism Inequalities).

Let $(H,Q,\Gamma)$ be a TAP instance with $H=(V,E)$ and
$\Gamma=G^{1},\ldots,G^{L}$. Further, let
$\tilde{G}\coloneqq\bigcup_{t=1}^{L}\tilde{G}^{t}$, $\tilde{G}^{t}\subseteq
G^{t}$, be a set of gates and let $C\coloneqq C(\tilde{G})=(\tilde{Q},E_{C})$
be the corresponding connectivity graph. Finally, let $d\geq 0$ be an integer
and let $H^{d}$ be the $d$-th relaxed graph of $H$. Now consider the time
steps $t_{0}\coloneqq\min\\{t\mid\tilde{G}^{t}\neq\emptyset\\}$ and
$t_{1}\coloneqq\max\\{t\mid\tilde{G}^{t}\neq\emptyset\\}$. If there is no
edge-induced subgraph of $H^{d}$ isomorphic to $C$, then the inequalities

$\displaystyle\sum_{t=t_{0}}^{t_{1}-1}\sum_{q\in\tilde{Q}}\sum_{(i,j)\in
V\times V}d_{H}(i,j)x^{t}_{q,i,j}$ $\displaystyle\geq d+1$ (7) and
$\displaystyle\sum_{t=t_{0}}^{t_{1}-1}\sum_{q\in Q}\sum_{(i,j)\in V\times
V}d_{H}(i,j)x^{t}_{q,i,j}$ $\displaystyle\geq 2(d+1)$ (8)

hold.

###### Proof.

First, we show the validity of (7). For a feasible binary solution to Model
(3.2) consider the allocation $a_{t_{0}}$ of qubits at time step $t_{0}$. From
$a_{t_{0}}$ construct the edge-induced subgraph
$H^{\prime}=(V^{\prime},E^{\prime})$ of $H^{d}=(V,E^{d})$ defined by the edge
subset

$\displaystyle E^{\prime}\coloneqq\left\\{\\{i,j\\}\in E^{d}\mid
a_{t_{0}}^{-1}(i),a_{t_{0}}^{-1}(j)\in\tilde{Q}\wedge\\{a_{t_{0}}^{-1}(i),a_{t_{0}}^{-1}(j)\\}\in
E_{C}\right\\}.$

The mapping defined by the allocation $a_{t_{0}}$ is not an isomorphism from
$H^{\prime}$ to $C$, since $C$ is not subgraph isomorphic to $H^{d}$ by
assumption. Furthermore, the definition of $E^{\prime}$ implies that for every
edge $\\{i,j\\}$ in $H^{\prime}$ there is a corresponding edge
$\\{a_{t_{0}}^{-1}(i),a_{t_{0}}^{-1}(j)\\}$ in $C$. This means that there is
an edge $\\{q_{1},q_{2}\\}$ in $C$ such that
$\\{a_{t_{0}}(q_{1}),a_{t_{0}}(q_{2})\\}$ is not an edge in $H^{\prime}$. This
edge corresponds to a gate $(q_{1},q_{2})\in G^{t^{\prime}}$ for some
$t^{\prime}$ with $t_{0}<t^{\prime}\leq t_{1}$. For this gate, the logical
qubits $q_{1}$ and $q_{2}$ are located at physical qubits more than $d+1$
apart at $t_{0}$, i.e.

$d_{H}(a_{t_{0}}(q_{1}),a_{t_{0}}(q_{2}))>d+1,$ (9)

but need to be at neighbouring physical qubits at $t^{\prime}$. Thus, together
they need to move a distance of at least $d+1$:

$d_{H}(a_{t_{0}}(q_{1}),a_{t^{\prime}}(q_{1}))+d_{H}(a_{t_{0}}(q_{2}),a_{t^{\prime}}(q_{2})\geq
d+1.$ (10)

This implies (7). Now, we show the validity of (8). For a subset of vertices
$\tilde{Q}\subseteq Q$ let

$L(\tilde{Q})\coloneqq\sum_{q\in\tilde{Q}}d_{H}\left(a_{t_{0}}(q),a_{t^{\prime}}(q)\right)$
(11)

denote the total distance moved by logical qubits in $\tilde{Q}$ between
$t_{0}$ and $t^{\prime}$. For the left-hand side of (7) it then holds

$\sum_{t=t_{0}}^{t_{1}-1}\sum_{q\in\tilde{Q}}\sum_{(i,j)\in V\times
V}d_{H}(i,j)x^{t}_{q,i,j}\geq L(\tilde{Q}),$ (12)

since $t^{\prime}\leq t_{1}$. Consider the $Q$-permutation $\pi\coloneqq
a_{t_{0}}^{-1}\circ a_{t^{\prime}}\colon Q\rightarrow Q$. Let
$K_{1},K_{2}\subseteq Q$ be the cycles in $\pi$ containing $q_{1}$ and
$q_{2}$, respectively. First, consider the case $K_{1}\neq K_{2}$. Here it
holds

$L(\tilde{Q})\geq L(K_{1})+L(K_{2})$ (13)

since $K_{1}\cap K_{2}=\emptyset$. Furthermore, we know $L(K_{i})\geq
2d_{H}(a_{t_{0}}(q_{i}),a_{t^{\prime}}(q_{i}))$, $i\in\\{1,2\\}$, since
$K_{1}$ and $K_{2}$ are cycles and $d_{H}(\cdot,\cdot)$ satisfies the triangle
inequality, see Figure 6a for an illustration.

(a) $K_{1}\neq K_{2}$

(b) $K_{1}=K_{2}$

Figure 6: Sketch for the two cases in the proof of Lemma 3.10. The cycle in
(a) contains exactly one of the qubits $q_{1}$ and $q_{2}$ and its length is
greater or equal than twice the distance this qubit moves (dashed line). The
cycle in (b) contains both qubits $q_{1}$ and $q_{2}$, and its length is
greater or equal than twice the initial distance of $q_{0}$ and $q_{1}$
(dashed line).

Thus, we have

$L(K_{1})+L(K_{2})\geq
2d_{H}(a_{t_{0}}(q_{1}),a_{t^{\prime}}(q_{1}))+2d_{H}(a_{t_{0}}(q_{2}),a_{t^{\prime}}(q_{2}))\overset{\eqref{eq:dist1q2traveled}}{\geq}2(d+1).$
(14)

On the other hand, if $K_{1}=K_{2}\eqqcolon K$ we find

$L(\tilde{Q})\geq L(K).$ (15)

Thus, again using the triangle inequality and the fact that $K$ is a cycle, we
have

$L(K)\geq
2d_{H}(a_{t_{0}}(q_{1}),a_{t_{0}}(q_{2}))\overset{\eqref{eq:distq1q2t0}}{\geq}2(d+1)$
(16)

as illustrated in Figure 6b. This implies (8). ∎

Clearly, any valid subgraph isomorphism inequality cuts off all solutions of
the LP relaxation with zero cost. However, generating these inequalities can
be computationally expensive, since subgraph isomorphism is NP-complete.

Having shown its NP-hardness and analysed its binary model, we conclude this
section about the TAP by introducing methods to reduce the model size for
larger instances.

###  3.4 Algorithmic Enhancements for Large Instances

To keep the size of Model (3.2) tractable for larger hardware graphs and
quantum circuits, we describe two heuristic methods for eliminating variables,
which means that both methods may in principle lead to suboptimal solutions.
However, our numerical tests show no or only marginal suboptimality together
with a significant reduction in computation time. The number of flow variables
is decreased by either removing edges from the time-expanded graph $G$ or by
removing logical qubits not participating in two-qubit gates.

##### Limiting Distance.

The first method eliminates edges from $G$. We only consider edges which
connect physical qubits in $H$ that are within a given distance limit. Thus,
we limit the distance that a logical qubit may move between subsequent
allocations.

##### Limiting the Number of Qubits.

In practice, the number of logical qubits is usually smaller than the number
of physical qubits. In principle, this can be reduced to the case of an equal
number of logical and physical qubits by adding inactive logical qubits not
participating in any two-qubit gates. However, the number of variables in the
TAP model grows linearly with $\lvert Q\rvert$. Therefore, we only consider
flow variables associated with logical qubits. This bears the problem that
costs associated with the flow of inactive qubits are not counted. To
circumvent this issue, we do not allow logical qubits to change positions with
inactive qubits. Thus, the subgraph of $H$ at which active logical qubits are
located is fixed for all time steps.

##  4 Solving the Token Swapping Problems

We now turn to the solution of the token swapping problems which need to be
solved as the second step in the proposed routing procedure once an optimal
TAP solution has been found. Indeed, we will describe two solution methods for
the token swapping problem: an efficient approximation algorithm as well as an
exact branch-and-bound approach.

###  4.1 Approximate Solution

The proposed qubit routing procedure employs a modified version of the
efficient token swapping algorithm given by Miltzow et al. in [16]. It has a
guaranteed approximation factor of four in terms of swap count. However, we
note that on all instances we considered the algorithm shows an approximation
factor better than 1.5 in practice. Indeed, the analysis in [16] shows that
the approximation factor is strictly less than four.

We briefly review the original approximation algorithm. For a detailed
description and analysis, we refer the interested reader to [16]. Given a
token swapping instance $(H,Q,a,a^{\prime})$, the algorithm iteratively
performs swaps by gradually changing the allocation from $a$ to $a^{\prime}$.
We adapt the notion of [16] for an arbitrary token allocation $\tilde{a}$: an
_unsatisfied vertex_ is a vertex $i\in H$ which does not hold its target
token, i.e. $\tilde{a}^{-1}(i)\neq(a^{\prime})^{-1}(i)$. A _distance-
decreasing neighbour_ of a vertex $i\in H$ is a neighbour of $i$ from which
the distance to the target vertex of the token sitting at $i$ is strictly less
than the distance from $i$. A _happy swap chain_ of length $k$ is a sequence
of $k$ swaps
$(i_{1},i_{2}),(i_{2},i_{3}),\ldots,(i_{k-1},i_{k}),(i_{k},i_{k+1})\subseteq\mathcal{T}_{H}$
such that every token involved has its distance decreased after the swaps in
the chain have been performed sequentially. An _unhappy swap_ is a swap that
moves a token already positioned at its target vertex. Given these notions,
the algorithm works as follows:

  1. 1.

Choose an unsatisfied vertex.

  2. 2.

Perform a walk along distance-decreasing neighbours until a cycle is closed or
a dead end encountered, i.e. no distance-decreasing neighbour exists.

  3. 3.

Perform the happy swap chain or unhappy swap.

  4. 4.

Go to 1.

The key insight why this algorithm terminates is that vertices involved in an
unhappy swap will be part of a future happy swap chain. Furthermore, the
number of happy swap chains is finite since they necessarily move tokens
further towards their targets. We modify the algorithm in two ways aiming to
further reduce both resulting swap count and depth. First, in step 1 we choose
an unsatisfied vertex that was not part of the previously performed swap chain
or unhappy swap. With this choice, we try to reduce depth: only swaps on
disjoint vertices can be performed in parallel. Second, in step 2 we try to
choose cleverly among all distance-decreasing neighbours by exploring them
first. If possible, we choose the one closing the smallest cylce. The
motivation for this is that all swaps in a happy swap chain cannot be
performed in parallel since they share common vertices. If no cycle can be
closed, we try to avoid dead ends, i.e. unhappy swaps. Intuitively, a typical
optimal solution does not use many unhappy swaps.

Having introduced an efficient but approximate method, we now derive an exact
branch and bound algorithm.

###  4.2 Exact Solution

The exact approach will be used for benchmarking the approximation algorithm
presented in the previous section. The exact method finds a shortest path from
the initial allocation to the final allocation in the associated _Cayleigh
graph_. The latter possesses a node for each of the $\lvert Q\rvert!$ possible
allocations. Two nodes are adjacent if and only if there is a swap, i.e. a
transposition $\tau\in\mathcal{T}_{H}$, which maps between the associated
allocations. It follows that an optimal solution to a token swapping problem
is given by a shortest path from the initial allocation to the final
allocation in the Cayleigh graph. The search algorithm works through the nodes
of the Cayleigh graph starting at the initial allocation. At each node, we
give upper and lower bounds on the shortest path length to the goal node.
Upper bounding is achieved by running the modified 4-approximation algorithm
from Section 4.1 yielding a path from the current node to the goal node. This
path is augmented by the path from the start node to the current node, giving
a feasible solution to the token swapping problem.

We use two different types of lower bounds, described in the following
sections.

####  4.2.1 Improved Distance Lower Bound

As already mentioned in Section 3, a well-known (cf. e.g. [28]) lower bound on
the optimal value of a token swapping instance $I=(H,Q,a,a^{\prime})$ is given
by

$\displaystyle\text{OPT(I)}\geq\sum_{q\in
Q}\frac{d_{H}(a(q),a^{\prime}(q))}{2}.$ (17)

Our aim is to improve upon this lower bound. To this end, we introduce the
notion of a _blocking vertex_.

######  Definition 4.1 (Blocking vertex).

Let $(H,Q,a,a^{\prime})$ be a token swapping instance. Let $q\in Q$ be a token
not yet sitting at its target vertex, i.e. $a(q)\neq a^{\prime}(q)$. A vertex
$v\in H$ is called _$q$ -blocking vertex_ if it is contained in a shortest
path in $H$ from $a(q)$ to $a^{\prime}(q)$, and if it holds its desired token
i.e. $a^{-1}(v)=(a^{\prime})^{-1}(v)$.

Intuitively, the presence of blocking vertices increases the number of
required swaps. This intuition is confirmed by

######  Lemma 4.2 (Improved distance lower bound).

Let $I=(H,Q,a,a^{\prime})$ be a token swapping instance. Let $Q^{*}\subseteq
Q$ denote the set of tokens not yet sitting at their target vertices. For
every $q\in Q^{*}$ let $B_{q}\subset V$ be a set of $q$-blocking vertices (see
Definition 4.1) of size $k_{q}$ such that

  1. (i)

any path in $H$ from $q$ to its target vertex avoiding any of the blocking
vertices in $B_{q}$ has at least length $d_{H}(a(q),a^{\prime}(q))+2k_{q}$ and

  2. (ii)

the sets $B_{q}$ are pairwise disjoint.

Then it holds

$\displaystyle\mathrm{OPT}(I)\geq\left\lceil\sum_{q\in
Q}\frac{d_{H}(a(q),a^{\prime}(q))}{2}\right\rceil+\sum_{q\in Q^{*}}k_{q}.$
(18)

Figure 7: Sketch for the proof of Lemma 4.2. The shortest path from the vertex
holding token $q$ to its target vertex (horizontal arrow) is blocked by
vertices already holding their desired tokens (black circles). In any
solution, $q$ takes a detour (curved arrow) or satisfied tokens at blocking
vertices are moved. Either alternative results in additional swaps.

###### Proof.

Consider a solution $\mathcal{S}=\tau_{1},\tau_{2},\ldots,\tau_{N}$ and denote
$[N]\coloneqq\\{1,\ldots,N\\}$. We partition
$Q^{*}=Q^{*}_{1}\,\dot{\cup}\,Q^{*}_{2}$ according to the following scheme:
for $q\in Q^{*}$ let $q\in Q^{*}_{1}$ if $q$ the path of $q$ in $\mathcal{S}$
contains all associated blocking vertices $B_{q}$, otherwise let $q\in
Q^{*}_{2}$.

For $q\in Q_{2}^{*}$ and $v,v^{\prime}\in V$ let
$\mathcal{P}_{q}(v,v^{\prime})$ denote the set of $v$-$v^{\prime}$-paths in
$H$ which do not contain all of the blocking vertices $B_{q}$. For a path
$p\in\mathcal{P}_{q}(v,v^{\prime})$ set $l(p)\coloneqq\lvert p\rvert-1$ and
define

$\displaystyle
d_{H,q}(v,v^{\prime})\coloneqq\min\\{l(p)\mid\in\mathcal{P}_{q}(v,v^{\prime})\\}$
(19)

as the edge distance between $v$ and $v^{\prime}$ if only paths not containing
all $q$-blocking vertices are allowed. Note that
$\mathcal{P}_{q}(a(q),a^{\prime}(q))\neq\emptyset$: from the definition of
$Q^{*}_{2}$ it follows that for all $q\in Q^{*}_{2}$ there is a path not
containing $B_{q}$. Furthermore, by property (i) of $B_{q}$ it holds

$d_{H,q}(a(q),a^{\prime}(q))\geq d_{H}(a(q),a^{\prime}(q))+2k_{q},$ (20)

for all $q\in Q^{*}_{2}$, see Figure 7 for an illustration.

For $i\in[N]$, consider the sum $d_{i}$ of all edge distances from tokens to
their targets after swap $\tau_{i}$ has been performed, accounting for the
detours taken by tokens in $Q^{*}_{2}$,

$\displaystyle d_{i}\coloneqq$ $\displaystyle\sum_{q\in
Q^{*}_{1}}d_{H}(\tau_{i}\circ\tau_{i-1}\circ\ldots\circ\tau_{1}\circ
a(q),a^{\prime}(q))+$ $\displaystyle\sum_{q\in
Q^{*}_{2}}d_{H,q}(\tau_{i}\circ\tau_{i-1}\circ\dots\circ\tau_{1}\circ
a(q),a^{\prime}(q)).$

Accordingly, set

$\displaystyle d_{0}\coloneqq\sum_{q\in
Q_{1}^{*}}d_{H}(a(q),a^{\prime}(q))+\sum_{q\in
Q^{*}_{2}}d_{H,q}(a(q),a^{\prime}(q)).$ (21)

Using (20) yields

$d_{0}\geq\sum_{q\in Q^{*}}d_{H}(a(q),a^{\prime}(q))+\sum_{q\in
Q^{*}_{2}}2k_{q}=\underbrace{\sum_{q\in
Q}d_{H}(a(q),a^{\prime}(q))}_{\eqqcolon D_{0}}+\sum_{q\in Q^{*}_{2}}2k_{q},$
(22)

where we used that $d_{H}(a(q),a^{\prime}(q))=0$ for all $q\in Q\setminus
Q^{*}$. Furthermore, for $i\in[N]$ let $\Delta d_{i}$ denote the difference

$\displaystyle\Delta d_{i}\coloneqq d_{i}-d_{i-1}.$ (23)

Clearly, $d_{N}=0$, since $\tau_{N}\circ\ldots\circ\tau_{2}\circ\tau_{1}\circ
a=a^{\prime}$. Thus,

$\sum_{i\in[N]}\Delta d_{i}=d_{N}-d_{0}=-d_{0}.$ (24)

Furthermore, for any $i\in[N]$ it holds $-2\leq\Delta d_{i}\leq 2$, since a
swap moves two tokens each for a distance of 1.

Now we turn to $Q_{1}^{*}$. By definition of $Q_{1}^{*}$, a vertex
$v\in\bigcup_{q\in Q^{*}_{1}}B_{q}$ is necessarily involved in a swap. Denote
$i_{v}\coloneqq\min\\{i\in[N]\mid v\in\tau_{i}\\}$ and let
$I\coloneqq\\{i_{v}\mid v\in\bigcup_{q\in Q^{*}_{1}}B_{q}\\}$. Partition
$I=I_{0}\,\dot{\cup}\,I_{2}$ in the following way. For $i_{v}\in I$, if there
is another $v^{\prime}\in\bigcup_{q\in Q^{*}_{1}}B_{q}$ such that
$i_{v}=i_{v^{\prime}}$ we let $i_{v}$ belong to the set $I_{2}$. Otherwise, we
let $i_{v}$ belong to the set $I_{0}$. Then, for $i\in I_{0}$ we have $\Delta
d_{i}\geq 0$, since at least one token leaves its target vertex by swap
$\tau_{i}$. Furthermore, for $i\in I_{2}$ we have $\Delta d_{i}=2$, since both
tokens leave their target vertices by swap $\tau_{i}$.

Let $n_{0}\coloneqq\lvert I_{0}\rvert$ and $n_{2}\coloneqq\lvert I_{2}\rvert$.
Using property (ii) of $B_{q}$ and the definition of $I_{0}$ and $I_{2}$, it
holds

$\displaystyle\lvert\bigcup_{q\in Q^{*}_{1}}B_{q}\rvert=\sum_{q\in
Q^{*}_{1}}k_{q}=2n_{2}+n_{0}.$ (25)

Let $\bar{I}\coloneqq[N]\setminus{I}=[N]\setminus(I_{0}\cup I_{2})$. Then we
can use (24) to obtain

$\displaystyle-d_{0}=\sum_{i\in[N]}\Delta d_{i}$
$\displaystyle=\underbrace{\sum_{i\in I_{0}}\Delta d_{i}}_{\geq
0}+\underbrace{\sum_{i\in I_{2}}\Delta
d_{i}}_{=2n_{2}}+\sum_{i\in\bar{I}}\Delta d_{i}\geq
2n_{2}+\sum_{i\in\bar{I}}\Delta d_{i}.$ (26)

Thus, with (22), we have

$\displaystyle\sum_{i\in\bar{I}}\Delta d_{i}\leq-D_{0}-\sum_{q\in
Q^{*}_{2}}2k_{q}-2n_{2}.$ (27)

Since $\Delta d_{i}\geq-2$, it follows

$\displaystyle\lvert\bar{I}\rvert\geq\left\lceil\frac{D_{0}+\sum_{q\in
Q^{*}_{2}}2k_{q}+2n_{2}}{2}\right\rceil=\left\lceil\frac{D_{0}}{2}\right\rceil+\sum_{q\in
Q^{*}_{2}}k_{q}+n_{2}.$ (28)

Altogether, this means

$\displaystyle N$ $\displaystyle=\lvert I_{0}\rvert+\lvert
I_{2}\rvert+\lvert\bar{I}\rvert\geq
n_{0}+n_{2}+\left\lceil\frac{D_{0}}{2}\right\rceil+\sum_{q\in
Q^{*}_{2}}k_{q}+n_{2}$
$\displaystyle=\left\lceil\frac{D_{0}}{2}\right\rceil+\sum_{q\in
Q^{*}_{2}}k_{q}+\sum_{q\in
Q^{*}_{1}}k_{q}=\left\lceil\frac{D_{0}}{2}\right\rceil+\sum_{q\in Q^{*}}k_{q}$

∎

As it turns out computationally, this lower bound is particularly strong on
graphs with low connectivity. Intuitively, a blocked path on a sparse graph
results in large detours. On the contrary, the second lower bound employed by
the branch-and-bound algorithm is stronger on dense graphs.

####  4.2.2 Complete Split Graph Lower Bound

While the bound from Lemma 4.2 is based on relaxing the restriction that
tokens need to be moved by swaps, we now derive a lower bound that is based on
relaxing the set of allowed transpositions $\mathcal{T}_{H}$. This set is
increased by adding edges to $H$. We use a result of Yasui et al. ([29]), who
give an efficient algorithm for solving token swapping problems on _complete
split graphs_.

######  Definition 4.3 (Complete split graph).

A vertex subset of a graph is called _independent_ if no two vertices in the
subset are connected by an edge. A vertex subset is called a _clique_ if any
two vertices in the subset are connected by an edge. A graph is called a
_complete split graph_ if it can be partitioned into an independent set and a
clique such that every vertex in the independent set is connected to all
vertices in the clique.

######  Lemma 4.4 (Complete split graph lower bound [29]).

Let $I=(H,Q,a,a^{\prime})$ be a token swapping instance. Furthermore, let
$V^{\prime}\subset V$ be an independent set in $H$. Let $n=\lvert V\rvert$ be
the number of vertices, let $r$ be the number of cycles in the permutation
$a\circ(a^{\prime})^{-1}$ including fix points, i.e. cycles of length 1. Let
$q$ be the number of cycles in the permutation $a\circ(a^{\prime})^{-1}$
excluding fix points that only involve vertices in $V^{\prime}$. Then it holds

$\displaystyle\mathrm{OPT}(I)\geq n-r+2q.$ (29)

###### Proof.

Consider the graph $H^{\prime}=(V,E^{\prime})$ which is derived from $H$ and
$V^{\prime}$ in the following way: connect every vertex in $V^{\prime}$ with
all vertices in $V\setminus V^{\prime}$. Additionally, add all missing edges
in $V\setminus V^{\prime}$. Then $H^{\prime}$ is a complete split graph.
Furthermore, $E\subseteq E^{\prime}$ holds, and thus also
$\mathcal{T}_{H}\subseteq\mathcal{T}_{H}^{\prime}$. Consider a solution
$\mathcal{S}$ of size $N$ to $I$. This is also a solution to the token
swapping instance
$I^{\prime}=(H^{\prime},Q,a\circ(a^{\prime})^{-1},\mathrm{id})$. Thus, we have

$\displaystyle\lvert\mathcal{S}\rvert\geq\mathrm{OPT}(I^{\prime}).$ (30)

According to [29], the optimal value for $I^{\prime}$ is given by

$\displaystyle\mathrm{OPT}(I^{\prime})=n-r+2q.$ (31)

Thus,

$\displaystyle\lvert\mathcal{S}\rvert\geq n-r+2q.$ (32)

∎

We remark that Lemma 4.4 extends the well-known fact that token swapping on
complete graphs has optimal value $n-r$ (cf. e.g. [13]).

As already mentioned, both lower bounds complement each other: while the
distance bound is strong on sparse graphs, the complete-split-graph bound is
strong on dense graphs. The following lemma allows to increase the lower
bounds from Lemmas 4.2 and 4.4 by one in some cases.

######  Lemma 4.5 (Parity of solutions).

Let $I=(H,Q,a,a^{\prime})$ be a token swapping instance. Let
$\mathrm{sgn}(\pi)$ denote the sign of a permutation $\pi$. For any solution
$\mathcal{S}=\tau_{1},\tau_{2},\ldots,\tau_{N}$ to $I$, it holds

$\displaystyle N\equiv\left[1-\mathrm{sgn}(a^{\prime}\circ
a^{-1})\right]/2\mod 2$ (33)

###### Proof.

For any solution $\mathcal{S}$, it holds

$\displaystyle a^{\prime}=\tau_{N}\circ\ldots\circ\tau_{2}\circ\tau_{1}\circ
a$ (34)

or, equivalently

$\displaystyle a^{\prime}\circ
a^{-1}=\tau_{N}\circ\ldots\circ\tau_{2}\circ\tau_{1}.$ (35)

Using the fact that
$\mathrm{sgn}(\sigma\circ\pi)=\mathrm{sgn}(\sigma)\cdot\mathrm{sgn}(\pi)$ for
any two permutations $\sigma,\pi$, we have

$\displaystyle\mathrm{sgn}(a^{\prime}\circ
a^{-1})=\mathrm{sgn}(\tau_{N})\cdots\mathrm{sgn}(\tau_{2})\cdot\mathrm{sgn}(\tau_{1})=(-1)^{N}.$

Thus, $N$ is either even or odd, independent of $S$, i.e.

$\displaystyle N\equiv\left[1-\mathrm{sgn}(a^{\prime}\circ
a^{-1})\right]/2\mod 2.$ (36)

∎

With these theoretical properties of our solution method to the qubit routing
by swap insertion problem established, we now present various experimental
results.

##  5 Experimental Results

In the following, we benchmark the proposed routing algorithm. We perform
experiments covering four aspects:

  1. 1.

The computation time and solution quality when solving the TAP model with
strengthening inequalities.

  2. 2.

The quality and speed of the token swapping approximation algorithm.

  3. 3.

The quality of the routing solution produced by the combined algorithm
compared to other routing methods.

  4. 4.

The quality of solutions produced by quantum algorithms compiled with our
routing method, executed on actual quantum hardware.

The datasets generated and analyzed during the current study are available
from the corresponding author on reasonable request.

###  5.1 Benchmarking the Strengthened TAP Model

In the first step of the proposed decomposition approach a TAP needs to be
solved. Here, we give numerical results showing how the introduction of
subgraph isomorphism cutting planes (Section 3.3) accelerates the solution of
the TAP binary linear program. The flow model (3.2) is implemented and solved
via the Python API of _Gurobi_ ([10]). Before parsing the model to the solver,
we generate a set of valid inequalities described in Section 3.3. These are
added to a cut pool by a callback routine that is invoked whenever a new
feasible solution is found. All inequalities satisfied with equality by the
new solution are added to the pool. The solver can add cuts from the pool to
the model in order to cut off relaxed solutions.

The set of valid inequalities is generated as follows: first, the graphs
$H^{d}$ for $0\leq d\leq\max_{i,j\in H}d_{H}(i,j)$ are created and transformed
into their respective line graphs. To generate a subset of gates for which a
valid cut may be derived, we first choose time steps $t_{0}$ and $t_{1}$ such
that $1\leq t_{0}<t_{1}\leq L$. Then, all gates in gate groups between $t_{0}$
and $t_{1}$ are considered. The corresponding connectivity graph is generated
and transformed into its line graph. To perform the node-induced-subgraph
isomorphism check, we employ an implementation of the VF2 algorithm (see [5,
4]) from the NetworkX python package ([11]). This procedure is repeated by
iterating over time steps $t_{0}$, $t_{1}$ until all possible time steps have
been checked or a time limit is reached.

In the first benchmark, we measure the computation time required for solving
TAP instances to optimality. We compare the total time taken with
strengthening inequalities to the time taken for the bare model (3.2). We
measure the ratio $r_{t}=t_{\mathrm{SGI}}/t_{\mathrm{bare}}$, where
$t_{\mathrm{SGI}}$ and $t_{\mathrm{bare}}$ are the runtimes of the
strengthened model and the bare model, respectively. The time limit for
cutting plane generation is set to ten seconds and included in the total
runtime $t_{\mathrm{SGI}}$. Based on practical experience, the value of ten
seconds is a good compromise between time investment and cut pool size.

In all instances $(H,Q,\Gamma)$, the hardware graph $H$ is an 8-vertex line
and $Q\coloneqq\\{q_{0}\ldots,q_{7}\\}$. Each of the gate groups in
$\Gamma=G^{1},G^{2},\ldots,G^{L}$ consists of 4 randomly chosen pairs of
qubits. This results in instances with as many gates per layer as possible,
making the instances computationally hard. The depth $L$ is varied from 7 to
10, where 10 instances are generated per value of $L$. The instance sizes
resemble the capabilities of currently available quantum hardware. All
instances are solved to optimality; the absolute solutions times lie between 6
s and 660 s. The results are given in Table 1a. First, we observe shorter
runtimes for the strengthened model on average for all instance sizes, with
savings of up to 46Ṡecond, we observe rather high standard deviations,
indicating a strong instance dependence.

For the same benchmark set, we measure the tightness of the lower bound on the
number of required swaps in an optimal routing solution, given by the optimal
TAP solution. To this end, the objective of the full routing model (BIP) from
[19] is adjusted to minimize swap count instead of depth. Solving this model
yields the optimal swap number. We measure the relative difference
$r_{b}=\frac{\mathrm{OPT}_{\mathrm{BIP}}-\mathrm{OPT}_{\mathrm{TAP}}}{\mathrm{OPT}_{\mathrm{BIP}}}$,
with $\mathrm{OPT}_{\mathrm{BIP}}$ and $\mathrm{OPT}_{\mathrm{TAP}}$ denoting
the optimal BIP solution and TAP solution, respectively. These results are
shown in Table 1a. We observe tight bounds with a relative difference of at
most 17% on average, together with a small standard deviation. This is a key
feature of our decomposition approach: optimal TAP solution values typically
give a good estimate for the number of required swaps.

$L$ | 7 | 8 | 9 | 10
---|---|---|---|---
$r_{t}$ | 0.77 $\pm$ 0.20 | 0.54 $\pm$ 0.18 | 0.62 $\pm$ 0.31 | 0.88 $\pm$ 0.17
$r_{b}$ | 0.14 $\pm$ 0.05 | 0.12 $\pm$ 0.06 | 0.17 $\pm$ 0.06 | 0.14 $\pm$ 0.02
(a) Runtime comparison and bound quality when solving to optimality.

$L$ | 11 | 12 | 13 | 14
---|---|---|---|---
$r_{g}$ | 0.76 $\pm$ 0.43 | 0.90 $\pm$ 0.59 | 0.85 $\pm$ 0.24 | 0.82 $\pm$ 0.14
(b) Optimality gap comparison when solving with limited runtime.

Table 1: Allocation benchmark results. The hardware graph is an 8-vertex line
and $\lvert Q\rvert=8$. 10 instances with 4 randomly chosen gates per layer
are generated per depth $L$. Each instance is solved with and without
additional strengthening inequalities. We compare the average ratio $r_{t}$ of
runtimes when solving the TAP instance to optimality in (a) and the average
ratio $r_{g}$ of optimality gaps after 300 s in (b). Additionally, for the
instances in (a) we give the relative difference $r_{b}$ between the optimal
TAP value and the number of required swaps in an optimal routing solution.
After $\pm$, we state the empirical standard deviation.

The second benchmark compares the optimality gap after a fixed runtime limit
of 300 s on larger instances of depths 11 to 14. Instances are generated the
same way as above. The _absolute optimality gap_ is defined as the difference
between the best solution value and the largest lower bound. The _relative
optimality gap_ is defined as the absolute optimality gap divided by the best
solution value. We measure the ratio $r_{g}$ between the relative optimality
gap of the strengthened model and the bare model. The results can be seen in
Table 1b. We observe a smaller gap for the strengthened model on average for
all instance sizes. An improvement of up to 24% is achieved.

To summarize, the separation of valid subgraph isomorphism inequalities
improves both runtime and optimality gap. Thus, also the algorithm for the
full routing problem benefits from generating valid inequalities for the TAP
subproblem.

###  5.2 Benchmarking the Token Swapping Approximation

The second step in the decomposition approach consists of several consecutive
token swapping instances. Here we evaluate the performance of the employed
token swapping approximation. All algorithms described in this section have
been implemented in Python from scratch. First, we compare the original
4-approximation algorithm given in [16] to our modified version. For this
purpose, we generate random token swapping instances, i.e. initial and final
allocations111Due to symmetry it is sufficient to only vary initial
allocations while fixing the final allocation. on three different types of
planar graphs with increasing _connectivity_. By connectivity we mean the
ratio of edges to vertices. As graph types we choose ring graphs, $2\times
n/2$-meshes (“ladder graphs”) and square meshes with node numbers of 16, 36
and 64 for each type. This choice of graph types and sizes aims at resembling
current realizations of quantum computer hardware connectivity graphs. We give
the average ratio $r_{\\#}$ of swap numbers between the modified version and
the original version as well as the average ratio of depths $r_{d}$ in Table
2.

| Ring | Ladder | Mesh
---|---|---|---
$N$ | $r_{\\#}$ | $r_{d}$ | $r_{\\#}$ | $r_{d}$ | $r_{\\#}$ | $r_{d}$
16 | 0.96 | 0.91 | 0.88 | 0.82 | 0.81 | 0.78
36 | 1.0 | 1.0 | 0.91 | 0.83 | 0.83 | 0.78
64 | 1.0 | 0.99 | 0.93 | 0.79 | 0.84 | 0.81
Table 2: Comparison of original and modified token swapping approximation
algorithms. We give the average ratio of inserted swaps
$r_{\\#}=\\#_{\text{mod}}/\\#_{\text{orig}}$ and resulting depths
$r_{d}=d_{\text{mod}}/d_{\text{orig}}$. Comparison are made on rings, ladders
and meshes with increasing vertex number $N$.

Firstly, note that these ratios are less than or equal to $1$ across all
instances. This means that our modification always performs as least as good
as the original version, or better (if the factor is smaller than $1$).
Furthermore, both swap and depth ratio are roughly constant across graph sizes
for a given graph type. However, the advantage of the modified version
increases with growing connectivity of the graph type. An intuitive
explanation of this behaviour is the increase of swapping possibilities with
growing edge number. This allows the modified version to employ its potential
of deciding for “good” swaps instead of choosing them randomly.

Secondly, we compare the modified 4-approximation algorithm to the exact
branch-and-bound algorithm in terms of runtime and solution quality. We
measure the approximation ratio $r_{\\#}=\\#_{\text{approx}}/\\#_{\text{opt}}$
as well as the runtime ratio $r_{t}=t_{\text{approx}}/t_{\text{bnb}}$ on
various instances. As graph types, we choose lines, rings and ladders. Again,
this choice tries to resemble real hardware graphs while we also want to make
use of the fact that token swapping can be solved efficiently on lines and
rings222In fact, token swapping can be solved efficiently on lines, rings,
cliques, stars and brooms, see [13].. This allows use to measure the
approximation quality for instances which are intractable for the exact
algorithm. Experimental results are summarized in Table 3.

First, we observe that for lines the approximation ratio is unity across all
instances. As it might be conjectured from these numerical results, one can
also convince oneself theoretically that the approximation algorithm always
finds an optimal solution on lines.

Furthermore, Table 3 shows that for ring and ladder graphs the average
approximation ratio is less than or equal to 1.25, which is well below the
theoretical value of 4. This indicates that the theoretical approximation
factor of 4 underestimates the approximation quality observed in practice. As
already mentioned, the proof in [16] delivers an approximation ratio of
strictly less than four.

Comparing runtimes, we note a drastic improvement with growing graph size on
all graph types. For the largest instances the exact algorithm was able to
solve in reasonable time, a reduction of up to 96.5% is achieved. Moreover,
all absolute runtimes for the approximation algorithm are well below 1 s.
Altogether, the approximation algorithm typically delivers close-to-optimal
solutions with a significant time saving on all instances considered.

| Line | Ring | Ladder
---|---|---|---
$N$ | $r_{\\#}$ | $r_{t}$ | $\bar{t}_{\text{appr}}$ | $r_{\\#}$ | $r_{t}$ | $\bar{t}_{\text{appr}}$ | $r_{\\#}$ | $r_{t}$ | $\bar{t}_{\text{appr}}$
5 | 1.0 | 0.2 | 0.0 | 1.25 | 0.3 | 0.0 | 1.04 | 0.6 | 0.0
10 | 1.0 | 0.05 | 0.0 | 1.19 | 0.2 | 0.0 | 1.19 | 0.2 | 0.0
12 | 1.0 | - | 0.0 | 1.23 | 0.08 | 0.0 | 1.13 | 0.035 | 0.0
50 | 1.0 | - | 0.8 | 1.17 | - | 0.4 |  |  | 0.4
Table 3: Comparison of token swapping approximation and exact algorithm. We
give the average approximation ratio
$r_{\\#}=\\#_{\text{appr}}/\\#_{\text{opt}}$ and runtime ratio
$r_{t}=t_{\text{appr}}/t_{\text{bnb}}$ as well as the average runtime of the
approximation algorithm $\bar{t}_{\text{appr}}$ in seconds. Instance graphs
are lines, rings and ladders of increasing vertex number $N$. Token swapping
can be solved efficiently on lines and rings, cf. [13]. This allows us to
compare approximation ratios even when the runtime is intractable for the
exact algorithm (indicated by “-”). For ladders, no efficient exact algorithm
exists. Therefore, no results are presented for large ladder graphs (empty
cells). Averages are taken over 100 random instances. Exceptions due to
excessive runtime are: 54 instances for $N=10$ lines, 22 for $N=12$ rings and
18 for $N=12$ ladders.

###  5.3 Benchmarking the Routing Algorithm

In the following, the proposed routing method (abbreviated Alloc+TS) is
compared to various routing algorithms integrated in different open source
quantum compilers: the binary program (BIP) by Nannicini et al. ([19]), Li et
al.’s Sabre heuristic ([15]), Qiskit’s default compiler Stochastic swap ([2]),
Cambridge Quantum Computing’s Tket ([6, 23]) as well as Google’s Cirq [8]. As
performance metrics, we measure the increase in gate number and depth relative
to the original circuit.

We briefly summarize the Python implementation of our routing algorithm:
first, the two-qubit gates of the input circuit are grouped in layers
containing as much gates as possible. These gate groups together with the
hardware graph form a TAP instance. Next, the network flow model of the TAP
instance is solved. The resulting TAP solution defines a set of token swapping
problems, which are then solved by the approximation algorithm described in
Section 4. Finally, the routed quantum circuit is constructed by changing the
target qubits in the gates of the original circuit according to the allocation
sequence and by adding the computed swaps between subsequent allocations. For
a given instance, we compare the circuit depths resulting from the different
routing procedures to the optimal depth, which is computed in advance.
Additionally, we compare the resulting gate counts. Benchmarks are performed
on two different classes of circuits. The first class is taken from a publicly
available benchmark library while the second is generated according to a
protocol used for benchmarking quantum hardware.

##### QUEKO circuits.

First, we use the publicly available benchmarks circuits generated by Tan and
Cong in [24]. They propose a method for generating random circuits on a given
hardware graph with given depth as well as given single and two-qubit gate
numbers. Their construction delivers instances of the routing by swap
insertion problem with optimal value $0$, i.e. there exists an allocation of
logical to physical qubits such that all gates can be performed without
inserting swaps.

We use the “benchmark for near-term feasibility” (BNTF) circuit set from [24].
It is subdivided into circuits on a 16-qubit hardware graph (Rigetti’s
Apsen-4) with few gates, and circuits on a 54 qubit graph (Google’s Sycamore)
with many gates. For each hardware graph, 90 circuits are available – 10 for
each depth in $\\{5,10,\ldots,45\\}$. On each instance the time limit for BIP
is the set to 3600 s. The results are visualized in Figures 8 (a)–(d). By
construction, our algorithm will find a zero-swap solution if the TAP is
solved to optimality. This can be observed in Figures 8a and 8c. For this
benchmark set, we achieve improvements in resulting depth of up to a factor of
2 on small circuits (Figure 8b) and up to factor of 5 on larger circuits
(Figure 8d), compared to the best heuristic. This improvement comes at the
cost of increased runtime. Indeed, the runtimes of all other heuristics
applied in this benchmark are in the order of a few seconds (not shown in
Figure 8), while the proposed algorithm takes up to 12 s on average for small
circuits (Figure 8b) and up to 5 min on average for larger circuits (Figure
8d). However, we note a drastic reduction of runtime compared to the exact
approach BIP (Figure 8b) which is unable to find optimal solutions on larger
instances in less than one hour.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 8: Comparison of different routing algorithms. Data points are averaged
over 10 instances. Error bars indicate the empirical standard deviation.
Instances are the BNTF circuits from [24] on Rigetti’s _Aspen-4_ 16-qubit
hardware graph (upper plots) as well as Google’s _Sycamore_ 54-qubit hardware
graph (middle plots). For every BNTF instance, a zero-swap solution exists.
Missing data points for BIP are due to excessive runtime. Instances for the
lower plots are QV circuits ([7]) of equal width and depth on an 8-qubit line
hardware graph. (a), (c), (e): Relative increase in the number of two-qubit
gates caused by swap insertion. (b), (d), (f): Relative depth increase as well
as runtime of Alloc+TS (black) and BIP (gray) on the second y-axis.

##### QV circuits.

For a second benchmark, we consider quantum volume (QV) circuits. According to
[17], a QV circuit of depth $d$ and width $m$ is a circuit on $m$ qubits with
$d$ layers where each layer consists of $\lfloor m/2\rfloor$ random two-qubit
gates. Thus, the resulting circuit is maximally dense in terms of gates per
layer. QV circuits of equal depth and width are used to benchmark quantum
computers (see [7]). We consider QV circuits of equal depth and width on an
8-qubit line hardware graph. We vary the depth from 4 to 8 and generate 10
instances per depth. Optimal-depth solutions are calculated using BIP. Our
results are visualized in Figures 8 (e) and (f). In Figure 8f, we observe that
the proposed algorithm finds solutions which are close-to-optimal in terms of
depth. Furthermore, it performs best among all considered routing heuristics.
This also holds when comparing gate numbers, as can be seen in Figure 8e.
Again, this improvement comes at the cost of increased runtime. While the
other heuristics require several seconds (not shown in Fig. 8), the proposed
method takes up to 10 minutes on average (Figure 8f). Still, the runtime is
significantly shorter than for the exact method BIP (Figure 8f).

We discuss three reasons why the improved routing solutions can be worth the
additional time investment. First, access to currently available quantum
hardware is often organized by a waiting queue, such that the time for
compilation is not a critical factor. Second, different quantum circuits with
equal connectivity graphs only need to be routed once. A popular example are
parameterized quantum circuits, which are often considered as promising
candidates for first quantum applications. The connectivity graph of such a
circuit is not parameter dependent. Third, current quantum hardware benefits
strongly from good routing solutions. This is demonstrated in the next
section.

###  5.4 Performance on Real Hardware

Lastly, we demonstrate the influence of routing quality on the performance of
current quantum hardware. For this purpose, we compile example quantum
algorithms to a real hardware device using different routing methods and
compare the results after execution.

The example algorithms are the Quantum approximate optimization algorithm
(QAOA) with depth parameter $p=1$ applied to MaxCut on ring graphs with 7, 8
and 9 vertices. QAOA is a well-known hybrid quantum-classical algorithm. A
detailed description is given by Farhi et al. in [9]. We choose QAOA for three
reasons: First, it is often considered as a promising candidate for near-term
quantum applications because of its universality and adaptive complexity.
Second, routing is an essential requirement when running QAOA on real
hardware, since QAOA employs two-qubit gates along the edges of an underlying
problem graph which usually is not isomorphic to the hardware connectivity
graph. This is also the case in our experiment. Third, the approximation ratio
of QAOA is an easily accessible, single-metric performance measure. The
approximation ratio is defined as $r=\mathbbm{E}(c)/c^{*}$. Here
$\mathbbm{E}(c)$ denotes the expected value of the cut size $c$ produced by
QAOA, and $c^{*}$ is the size of an optimal cut. QAOA is a parameterized
quantum algorithm. For QAOA with depth $p=1$ applied to MaxCut on ring graphs
the optimal parameters can be calculated analytically, cf. [25].

The routing algorithms tested are the proposed method Alloc+TS, Tket and
Stochastic swap. For circuit construction, compilation and backend
communication, we use IBMs SDK Qiskit ([2]). The quantum backend in our
experiment is _ibmq_ehnigen_ ([1]). We compile each of the QAOA circuits with
all three routing methods while other compiler options remain unchanged.

Table 4 compares the gate counts and depths of the compiled circuits as well
as the achieved approximation ratios after execution on _ibmq_ehnigen_. The
proposed routing procedure All+TS performs best both in terms of swaps added
and depth increase on all instances, followed by Tket. We observe that small
gate count and depth correlate with an increase of approximation ratio on all
instances. For the instance with $N=8$, we note large reductions in gate count
of 65% and 69% when compared to Tket and Stochastic swap, respectively.
Similarly, the circuit depth reductions amount to 60% and 68%, respectively.
For the measured approximation ratio, we observe an increase of 23% and 25%,
respectively. This example shows a clear correlation between routing quality
and quantum algorithm performance.

Although this experiment is by no means a comprehensive benchmark of our
routing method, it stresses the need for good compiling methods in quantum
computation on currently available hardware and shows clear benefits of our
approach.

$N$ | $\\#_{2q}^{\text{ATS}}$ | $\\#_{2q}^{\text{Tket}}$ | $\\#_{2q}^{\text{Sto}}$ | $d^{\text{ATS}}$ | $d^{\text{Tket}}$ | $d^{\text{Sto}}$ | $r_{\text{ideal}}$ | $r_{\text{ATS}}$ | $r_{\text{Tket}}$ | $r_{\text{St}}$
---|---|---|---|---|---|---|---|---|---|---
7 | 21 | 31 | 54 | 22 | 28 | 41 | 5.25 | 4.60 | 3.65 | 3.37
8 | 25 | 71 | 80 | 21 | 53 | 65 | 6.0 | 5.09 | 4.13 | 4.08
9 | 30 | 54 | 63 | 21 | 39 | 45 | 6.75 | 5.22 | 4.62 | 4.40
Table 4: Results for compiling quantum algorithms to real hardware with
different routing algorithms. Compiled quantum circuits are QAOA for MaxCut on
ring graphs with $N$ vertices. The target quantum device is _ibmq_ehningen_.
Given are the total two-qubit gate counts $\\#_{2q}^{\text{ATS}}$,
$\\#^{\text{Tket}}_{2q}$, $\\#_{2q}^{\text{Sto}}$ for the circuit routed with
Alloc+TS, Tket and Stochastic swap, respectively, as well as the corresponding
depths $d^{\text{ATS}}$, $d^{\text{Tket}}$, $d^{\text{Sto}}$. Furthermore, we
give the theoretical approximation ratio $r_{\text{ideal}}$ as well as the
actually achieved approximation ratios $r_{\text{ATS}}$, $r_{\text{Tket}}$,
$r_{\text{St}}$.

##  6 Conclusion

In this article, we have studied the problem of qubit routing by swap
insertion. This problem needs to be solved in the context of quantum
compilation but is practically intractable even for moderate instance sizes
when solved by an exact approach. We have proposed a new decomposition
approach that exploits the capabilities of exact integer programming while
reducing complexity significantly, at the expense of possibly returning
suboptimal solutions. Improvement in performance is achieved by solving only a
subproblem by integer programming. The objective of this subproblem, called
token allocation problem (TAP), is a lower bound on the overall routing
problem. We strengthened the model of the TAP by deriving valid inequalities.
The solution of the TAP defines the second subproblem, which is the well-known
problem of token swapping. Its solution yields the actual objective value of
the overall routing problem. For the token swapping problem, we have enhanced
an existing approximation algorithm and have developed an exact branch-and-
bound algorithm. Combining the solutions of both subproblems produces a
solution to the overall routing problem.

Finally, we have given experimental results showing the applicability of the
proposed algorithm. Although it being a heuristic approach, we observed close-
to-optimal solutions outperforming other heuristics. This improvement comes at
the expense of increased solution time; however, the solution time is still
much smaller than for exact approaches. Moreover, the lower bounds computed in
the TAP turn out to be tight, which we consider as a reason for the good
performance of the proposed decomposition method.

Our experiments on actual quantum hardware show the impact of high-quality
solutions to qubit routing for quantum computation on current hardware. For
the example algorithm QAOA, we observed an increase in approximation ratio
when compiled with the proposed routing method compared to standard
heuristics.

A possible field of future research is the structural analysis of the TAP in
order to further strengthen its integer programming formulation. Additionally,
practical experience with current quantum hardware suggests more complicated
cost functions to be considered: specific qubits and gates suffer from higher
error rates than others. Even more complicated, executing gates on specific
qubit pairs simultaneously results in high error rates – a phenomenon known as
_crosstalk_. Incorporating these effects in the routing algorithm is an
additional direction of future research.

## Acknowledgements

This research has been supported by the Bavarian Ministry of Economic Affairs,
Regional Development and Energy with funds from the Hightech Agenda Bayern as
well as the Center for Analytics – Data – Applications (ADA-Center) within the
framework of “BAYERN DIGITAL II” (20-3410-2-9-8).

## References

  * [1] IBM Quantum.  https://quantum-computing.ibm.com/, 2021.
  * [2] Qiskit: An open-source framework for quantum computing.  http://www.qiskit.org, 2021.  doi:10.5281/zenodo.2573505.
  * [3] Andrew M. Childs, Eddie Schoute, and Cem M. Unsal.  Circuit Transformations for Quantum Architectures.  In Wim van Dam and Laura Mancinska, editors, 14th Conference on the Theory of Quantum Computation, Communication and Cryptography (TQC 2019), volume 135 of Leibniz International Proceedings in Informatics (LIPIcs), pages 3:1–3:24, Dagstuhl, Germany, 2019. Schloss Dagstuhl–Leibniz-Zentrum für Informatik.  URL: http://drops.dagstuhl.de/opus/volltexte/2019/10395, doi:10.4230/LIPIcs.TQC.2019.3.
  * [4] Luigi P. Cordella, Pasquale Foggia, Carlo Sansone, and Mario Vento.  An improved algorithm for matching large graphs.  In 3rd IAPR-TC15 Workshop on Graph-based Representations in Pattern Recognition, Cuen, pages 149–159, 2001.
  * [5] Luigi P. Cordella, Pasquale Foggia, Carlo Sansone, and Mario Vento.  A (sub)graph isomorphism algorithm for matching large graphs.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 26(10):1367–1372, 2004.  doi:10.1109/TPAMI.2004.75.
  * [6] Alexander Cowtan, Silas Dilkes, Ross Duncan, Alexandre Krajenbrink, Will Simmons, and Seyon Sivarajah.  On the Qubit Routing Problem.  In Wim van Dam and Laura Mancinska, editors, 14th Conference on the Theory of Quantum Computation, Communication and Cryptography (TQC 2019), volume 135 of Leibniz International Proceedings in Informatics (LIPIcs), pages 5:1–5:32, Dagstuhl, Germany, 2019. Schloss Dagstuhl–Leibniz-Zentrum für Informatik.  URL: http://drops.dagstuhl.de/opus/volltexte/2019/10397, doi:10.4230/LIPIcs.TQC.2019.5.
  * [7] Andrew W. Cross, Lev S. Bishop, Sarah Sheldon, Paul D. Nation, and Jay M. Gambetta.  Validating quantum computers using randomized model circuits.  Physical Revue A, 100:032328, 2019.  URL: https://link.aps.org/doi/10.1103/PhysRevA.100.032328, doi:10.1103/PhysRevA.100.032328.
  * [8] Cirq Developers.  Cirq, August 2021.  See full list of authors on Github: https://github.com/quantumlib/Cirq/graphs/contributors.  doi:10.5281/zenodo.5182845.
  * [9] Edward Farhi, Jeffrey Goldstone, and Sam Gutmann.  A quantum approximate optimization algorithm, 2014.  arXiv:1411.4028.
  * [10] Gurobi Optimization, LLC.  Gurobi Optimizer Reference Manual, 2022.  URL: https://www.gurobi.com.
  * [11] Aric A. Hagberg, Daniel A. Schult, and Pieter J. Swart.  Exploring network structure, dynamics, and function using networkx.  In Gaël Varoquaux, Travis Vaught, and Jarrod Millman, editors, Proceedings of the 7th Python in Science Conference, pages 11 – 15, Pasadena, CA USA, 2008.  URL: http://conference.scipy.org/proceedings/SciPy2008/paper_2/.
  * [12] Scott Johnstun and Jean-François Van Huele.  Understanding and compensating for noise on IBM quantum computers.  American Journal of Physics, 89(10):935–942, 2021.  arXiv:https://doi.org/10.1119/10.0006204, doi:10.1119/10.0006204.
  * [13] Dohan Kim.  Sorting on graphs by adjacent swaps using permutation groups.  Computer Science Review, 22:89–105, 2016.  URL: https://www.sciencedirect.com/science/article/pii/S1574013716300077, doi:https://doi.org/10.1016/j.cosrev.2016.09.003.
  * [14] Aleks Kissinger and Arianne Meijer-Van De Griend.  CNOT circuit extraction for topologically-constrained quantum memories.  Quantum Information and Computation, 20(7/8):581–596, 2020.
  * [15] Gushu Li, Yufei Ding, and Yuan Xie.  Tackling the qubit mapping problem for NISQ-era quantum devices.  In Proceedings of the Twenty-Fourth International Conference on Architectural Support for Programming Languages and Operating Systems, ASPLOS ’19, page 1001–1014, New York, NY, USA, 2019. Association for Computing Machinery.  doi:10.1145/3297858.3304023.
  * [16] Tillmann Miltzow, Lothar Narins, Yoshio Okamoto, Günter Rote, Antonis Thomas, and Takeaki Uno.  Approximation and Hardness of Token Swapping.  In Piotr Sankowski and Christos Zaroliagis, editors, 24th Annual European Symposium on Algorithms (ESA 2016), volume 57 of Leibniz International Proceedings in Informatics (LIPIcs), pages 66:1–66:15, Dagstuhl, Germany, 2016. Schloss Dagstuhl–Leibniz-Zentrum für Informatik.  URL: http://drops.dagstuhl.de/opus/volltexte/2016/6408, doi:10.4230/LIPIcs.ESA.2016.66.
  * [17] Nikolaj Moll, Panagiotis Barkoutsos, Lev S. Bishop, Jerry M. Chow, Andrew Cross, Daniel J. Egger, Stefan Filipp, Andreas Fuhrer, Jay M. Gambetta, Marc Ganzhorn, Abhinav Kandala, Antonio Mezzacapo, Peter Müller, Walter Riess, Gian Salis, John Smolin, Ivano Tavernelli, and Kristan Temme.  Quantum optimization using variational algorithms on near-term quantum devices.  Quantum Science and Technology, 3(3):030503, 2018.  doi:10.1088/2058-9565/aab822.
  * [18] Lorenzo Moro, Matteo G. A. Paris, Marcello Restelli, and Enrico Prati.  Quantum compiling by deep reinforcement learning.  Communications Physics, 4(1):178, Aug 2021.  doi:10.1038/s42005-021-00684-3.
  * [19] Giacomo Nannicini, Lev S. Bishop, Oktay Gunluk, and Petar Jurcevic.  Optimal qubit assignment and routing via integer programming, 2021.  arXiv:2106.06446.
  * [20] Michael A. Nielsen and Isaac L. Chuang.  Quantum Computation and Quantum Information: 10th Anniversary Edition.  Cambridge University Press, 2010.  doi:10.1017/CBO9780511976667.
  * [21] Marcos Yukio Siraichi, Vinícius Fernandes dos Santos, Caroline Collange, and Fernando Magno Quintão Pereira.  Qubit allocation as a combination of subgraph isomorphism and token swapping.  Proc. ACM Program. Lang., 3(OOPSLA), Oct 2019.  doi:10.1145/3360546.
  * [22] Marcos Yukio Siraichi, Vinicius Fernandes dos Santos, Caroline Collange, and Fernando Magno Quintão Pereira.  Qubit Allocation.  In CGO 2018 - International Symposium on Code Generation and Optimization, pages 1–12, Vienna, Austria, Feb 2018.  URL: https://hal.archives-ouvertes.fr/hal-01655951, doi:10.1145/3168822.
  * [23] Seyon Sivarajah, Silas Dilkes, Alexander Cowtan, Will Simmons, Alec Edgington, and Ross Duncan.  t$|$ket$\rangle$: a retargetable compiler for NISQ devices.  Quantum Science and Technology, 6(1):014003, 2020.  doi:10.1088/2058-9565/ab8e92.
  * [24] Bochen Tan and Jason Cong.  Optimality study of existing quantum computing layout synthesis tools.  IEEE Transactions on Computers, 70(9):1363–1373, 2021.  doi:10.1109/TC.2020.3009140.
  * [25] Zhihui Wang, Stuart Hadfield, Zhang Jiang, and Eleanor Rieffel.  The quantum approximation optimization algorithm for MaxCut: A fermionic view.  Physical Review A, 97:022304, 2017.  doi:10.1103/PhysRevA.97.022304.
  * [26] Ingo Wegener.  Complexity Theory: Exploring the Limits of Efficient Algorithms, page 81.  Springer Berlin Heidelberg, 2005.  doi:https://doi.org/10.1007/3-540-27477-4.
  * [27] Robert Wille, Lukas Burgholzer, and Alwin Zulehner.  Mapping quantum circuits to IBM QX architectures using the minimal number of SWAP and H operations.  In Design Automation Conference, 2019.
  * [28] Katsuhisa Yamanaka, Erik D. Demaine, Takehiro Ito, Jun Kawahara, Masashi Kiyomi, Yoshio Okamoto, Toshiki Saitoh, Akira Suzuki, Kei Uchizawa, and Takeaki Uno.  Swapping labeled tokens on graphs.  Theoretical Computer Science, 586:81–94, 2015.  URL: https://www.sciencedirect.com/science/article/pii/S0304397515001656, doi:https://doi.org/10.1016/j.tcs.2015.01.052.
  * [29] Gaku Yasui, Kouta Abe, Katsuhisa Yamanaka, and Takashi Hirayama.  Swapping labeled tokens on complete split graphs.  IPSJ SIG Technical Report, 14, 2015.
  * [30] Alwin Zulehner, Alexandru Paler, and Robert Wille.  An efficient methodology for mapping quantum circuits to the ibm qx architectures.  IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems, 38(7):1226–1236, 2019.  doi:10.1109/TCAD.2018.2846658.


FTUV-2012-08-28 UAB-FT-719 IFIC-12-61

# MC generator TAUOLA: implementation of Resonance Chiral Theory for two and
three meson modes. Comparison with experiment

O. Shekhovtsova  I. M. Nugent  T. Przedzinski  P. Roig  Z. Was

###### Abstract

We present a partial upgrade of the Monte Carlo event generator TAUOLA with
the two and three hadron decay modes using the theoretical models based on
Resonance Chiral Theory. These modes account for 88% of total hadronic width
of the tau meson. First results of the model parameters have been obtained
using BaBar data for 3$\pi$ mode.

###### Keywords:

Tau physics, Monte Carlo event generator, chiral models

###### :

13.35.Dx; 12.39.Fe; 29.85.-c; 89.20.Ff

##  1 Introduction

Tau lepton is the only lepton that decays into hadrons. From the perspective
of low energies, as explored at B-factories, BaBar babar  and Belle belle
experiments, hadronic decay modes of tau lepton provide an excellent
laboratory for modeling hadronic interactions at the energy scale of about 1-2
GeV, where neither perturbative QCD methods nor chiral lagrangians are
expected to work with a good precision Braaten:1990ef . From the perspective
of high energy experiments such as those at LHC, good understanding of tau
leptons properties contributes important ingredients of new physics
signatures. With the discovery of a new particle around the mass of 125-126GeV
higgs:2012 , tau decays are an important ingredient for determining if this is
the Standard Model Higgs. This is especially pertinent since CMS has reported
a deficit in the number of fermion decays from the new particle relative to
the Standard Model Higgs Prediction.

The first version of the program TAUOLA Jadach:1993hs  was written in the
90’s. TAUOLA is Monte Carlo event generator which simulates tau decays for
both the leptonic and hadronic decays modes. The hadronic currents implemented
in TAUOLA are based on resonance dominance model (RDM) Kuhn  111Based on
Jadach:1993hs  Cleo and Aleph collaborations developed their own versions, non
published, which differ from one another by the intermediate resonance states
because of different detector sensitivity.. In the framework of RDM the
hadronic current for a three-pseudoscalar decay mode is a sum of weighted
products of the Breit-Wigner functions. It is demonstrated in GomezDumm:2003ku
that this approach is able to reproduce only LO $\chi$PT result. Moreover, the
model was not sufficient to describe the Cleo $KK\pi$ data Coan:2004ep . This
resulted in the Cleo collaboration reshaping the model by introducing two ad-
hoc parameters that spoilt the QCD normalization for Weiss-Zumino
contribution. The parameters are obtained fitting to data. However, before
making conclusion that the Wess-Zumino anomaly normalization is spoilt it
should be checked whether on oversimplified theoretical approximation, like
RDM, was applied.

The alternative approach based on the Resonance Chiral Theory (R$\chi$T) rcht
was proposed in GomezDumm:2003ku . The computations done within R$\chi$T are
able to reproduce the low-energy limit of $\chi$PT at least up to NLO and
demonstrate the right falloff in the high energy region Roig:talk . The
hadronic currents for the two and three meson modes ($\pi\pi$, $K\pi$, $KK$,
$\pi\pi\pi$ and $KK\pi$) calculated in the framework of R$\chi$T were
implemented in TAUOLA Shekhovtsova:2012ra . The tar-ball version of the code
is available at web-page web:RChL .

##  2 Results for two and three meson final states

It is of utmost importance to implement the theory in a form that is as useful
for general applications as possible. Therefore the corresponding hadronic
currents ($J^{\mu}$) in the upgraded version of TAUOLA have been written in
the most general form compatible with Lorentz invariance. For $\tau$ decay
channels with two mesons ($h_{1}(p_{1})$ and $h_{2}(p_{2})$)

$J^{\mu}=N\bigl{[}(p_{1}-p_{2}-\frac{\Delta_{12}}{s}(p_{1}+p_{2}))^{\mu}F^{V}(s)+\frac{\Delta_{12}}{s}((p_{1}+p_{2})^{\mu}F^{S}(s)\bigr{]},$
(1)

where $\Delta_{12}=m_{1}^{2}-m_{2}^{2}$, $s=(p_{1}+p_{2})^{2}$. For the final
state of three pseudoscalars, with momenta $p_{1}$, $p_{2}$ and $p_{3}$,
Lorentz invariance determines the decomposition

$\displaystyle J^{\mu}$ $\displaystyle=N$
$\displaystyle\bigl{\\{}T^{\mu}_{\nu}\bigl{[}c_{1}(p_{2}-p_{3})^{\nu}F_{1}+c_{2}(p_{3}-p_{1})^{\nu}F_{2}+c_{3}(p_{1}-p_{2})^{\nu}F_{3}\bigr{]}$
(2) $\displaystyle+c_{4}q^{\mu}F_{4}-{i\over
4\pi^{2}F^{2}}c_{5}\epsilon^{\mu}_{.\
\nu\rho\sigma}p_{1}^{\nu}p_{2}^{\rho}p_{3}^{\sigma}F_{5}\bigr{\\}},$

where: $T_{\mu\nu}=g_{\mu\nu}-q_{\mu}q_{\nu}/q^{2}$ denotes the transverse
projector, $q^{\mu}=(p_{1}+p_{2}+p_{3})^{\mu}$ is the momentum of the hadronic
system, $F$ is the pion decay constant in the chiral limit.

The decomposition (1) and (2) is the most general one, model-dependence is
included in the hadronic form factors ($F_{V}$, $F_{S}$ as well as $F_{i}$, i
= 1…5). The hadronic form factors calculated within R$\chi$T can be written as

$F_{I}=F_{I}^{\chi}+F_{I}^{R}+F_{I}^{RR}$ (3)

where $F_{I}^{\chi}$ is the chiral contribution, $F_{I}^{R}$ is the one
resonance contribution and $F_{I}^{RR}$ is the double-resonance part. The
explict form of the functions $F_{i}$ for 3$\pi$ and $KK\pi$ modes can be
found in Shekhovtsova:2012ra , Section 2, as well as in GomezDumm:2003ku ;
Dumm . The theoretical results for the hadronic currents were obtained in the
isospin limit ($m_{\pi}=(2m_{\pi^{+}}+m_{\pi^{0}})/3$,
$m_{K}=(m_{K^{+}}+m_{K^{0}})/2$), except for the two pion and two kaons modes.
In the phase space generation, the differences between neutral and charged
pion and kaon masses is taken into account, i.e. physical values are chosen.
This has to be done to obtain proper kinematic configurations. The model
parameters, more specifically the masses of the resonances and the coupling
constants, were fitted to Aleph data, requiring correct high-energy behaviour
of the related form factors, see Appendix C in Shekhovtsova:2012ra .

To check stability of multidimensional integration in TAUOLA the MC results
were compared with the semi-analytical ones (Gauss integration of analytical
results was used). The difference between MC prediction and semi-analytical
results for the partial decay width is less than $0.02\%$ for all channels.
Both differential spectra and numerical tests are collected at web:RChL .

In Table 1, the partial decay width values from PDG Nakamura:2010zzi  are
compared with our results obtained with isospin-averaged pseudoscalar masses
and with the physical ones. Comparison of the last two columns illustrates
numerical effect of physical masses. The difference between the MC result and
PDG one is $2\%-24\%$ depending on the channel. As expected, that agreement is
not good because only minimal attempts on adjusting to the model parameters
have been applied for the comparison with BaBar and Belle data. The next
section presents the first step toward this direction.

Channel | Width, [GeV]
---|---
| PDG | Equal masses | Phase space
|  |  | with masses
$\pi^{-}\pi^{0}\;\;\;\;$ | ($5.778\pm 0.35\%)\cdot 10^{-13}$ | ($5.2283\pm 0.005\%)\cdot 10^{-13}$ | $(5.2441\pm 0.005\%)\cdot 10^{-13}$
$K^{-}\pi^{0}\;\;\;\;$ | ($9.72\;\pm 3.5\%\;)\cdot 10^{-15}$ | ($8.3981\pm 0.005\%)\cdot 10^{-15}$ | $(8.5810\pm 0.005\%)\cdot 10^{-15}$
$\pi^{-}\bar{K}^{0}\;\;\;\;$ | ($1.9\;\;\;\pm 5\%\;\;\;)\cdot 10^{-14}$ | ($1.6798\pm 0.006\%)\cdot 10^{-14}$ | $(1.6512\pm 0.006\%)\cdot 10^{-14}$
$K^{-}K^{0}\;\;\;\;$ | ($3.60\;\pm 10\%\;\;)\cdot 10^{-15}$ | ($2.6502\pm 0.007\%)\cdot 10^{-15}$ | $(2.6502\pm 0.008\%)\cdot 10^{-15}$
$\pi^{-}\pi^{-}\pi^{+}$ | ($2.11\;\pm 0.8\%\;\;)\cdot 10^{-13}$ | ($2.1013\pm 0.016\%)\cdot 10^{-13}$ | $(2.0800\pm 0.017\%)\cdot 10^{-13}$
$\pi^{0}\pi^{0}\pi^{-}$ | ($2.10\;\pm 1.2\%\;\;)\cdot 10^{-13}$ | ($2.1013\pm 0.016\%)\cdot 10^{-13}$ | $(2.1256\pm 0.017\%)\cdot 10^{-13}$
$K^{-}\pi^{-}K^{+}$ | ($3.17\;\pm 4\%\;\;\;)\cdot 10^{-15}$ | ($3.7379\pm 0.024\%)\cdot 10^{-15}$ | $(3.8460\pm 0.024\%)\cdot 10^{-15}$
$K^{0}\pi^{-}\bar{K^{0}}$ | ($3.9\;\;\pm 24\%\;\;)\cdot 10^{-15}$ | ($3.7385\pm 0.024\%)\cdot 10^{-15}$ | $(3.5917\pm 0.024\%)\cdot 10^{-15}$
$K^{-}\pi^{0}K^{0}$ | ($3.60\;\pm 12.6\%\;\;)\cdot 10^{-15}$ | ($2.7367\pm 0.025\%)\cdot 10^{-15}$ | $(2.7711\pm 0.024\%)\cdot 10^{-15}$
Table 1: The $\tau$ decay partial widths. The PDG value Nakamura:2010zzi  (2nd
collumn) is compared with numerical results from TAUOLA with the R$\chi$T
currents: in the isospin limit for pseudoscalar masses (3rd collumn), using
physical masses (4th collumn).

##  3 Fit for three pion mode to BaBar data

The main problem with upgrading the MC simulation of $\tau$ decays is a lack
of the published spectrum. Currently, only the two pion modes belle  and three
pseudoscalar modes Nugent:2009zz  are published.

The results from a preliminary fit to the three pion mode
($\pi^{+}\pi^{-}\pi^{-}$) Nugent:2009zz  can be seen in Fig. 1 and Table 2.
Both the three particle and the $\pi^{+}\pi^{-}$ invariant mass distributions
have been considered. Disagreement about 12% level is visible in the low
energy region of two particle invariant mass distribution whereas for the
three invariant mass spectrum the difference between MC and data is less than
7%. However, the R$\chi$T parametrization is in better agreement with BaBar
data than CLEO one, see Fig. 1.

| $M_{\rho^{\prime}}$ | $\Gamma_{\rho^{\prime}}$ | $M_{a_{1}}$ | $F$ | $F_{V}$ | $F_{A}$ | $\beta_{\rho^{\prime}}$
---|---|---|---|---|---|---|---
Min. | 1.44 | 0.32 | 1.00 | 0.0920 | 0.12 | 0.1 | -0.36
Max | 1.48 | 0.39 | 1.24 | 0.0924 | 0.24 | 0.2 | -0.18
Default | 1.453 | 0.4 | 1.12 | 0.0924 | 0.18 | 0.149 | -0.25
Fit 222 | 1.4302 | 0.376061 | 1.21706 | 0.092318 | 0.121938 | 0.11291 | -0.208811
Table 2: Numerical values of the R$\chi$T parameters fitted to BaBar data for
three pion mode

The fit was done taking into account only $P$-wave contribution of two pion
system. As suggesting in Shibata:2002uv  the discrepance in the low mass
region could be described using a contribution from the scalar particle,
$S$-wave contribution. We expect that inclusion of $S$-wave contribution
Roig:talk  will improve the value of $\chi^{2}$. The values in the fifth row
of Table 2 are only the preliminary results. They will not necessarily
correspond to the minimum of $\chi^{2}/ndf$ of the final fit. Work is in
progress.

---
Figure 1: Invariant $\pi^{+}\pi^{-}\pi^{-}$ (up) and $\pi^{+}\pi^{-}$ (down)
mass distributions: the plots on the left-hand side correspond to the
differential decay distribution, the ones on the right-hand side to plot
ratios between MC and BaBar data. Lighter grey histograms are for R$\chi$T
parametrization, darker grey is for CLEO one.

##  4 Conclusion

The theoretical results for the hadronic currents of two and three
pseudoscalar modes, namely, $\pi\pi$, $K\pi$, $KK$, $\pi\pi\pi$ and $KK\pi$,
in the framework of R$\chi$T have been implemented in TAUOLA. These modes,
together with the one-meson decay modes, represent more than 88% of the
hadronic width of tau. R$\chi$T is a more controlled QCD-based model than the
usual used Breit-Wigner parametrizations. However, before making conclusion
about validity of the model the theoretical results have to be confronted with
the experimental data, which requires fit of the model parameters. Now that
the technical work on current installation is complete, the work on fits is in
progress in collaboration with theoreticians and experimentalists.

At LHC tau decays are only used for identification and are not used to study
their dynamic. However, the dynamics of tau decays are important for both
modeling the decays and -therefore the reconstruction and identification- and
for measuring the polarization of tau decays. Therefore, an upgrade to the
TAUOLA based on the BaBar and Belle results on tau decays is urgently needed
for LHC.

This research was supported by a Marie Curie Intra European Fellowship within
the 7th European Community Framework Programme (O.S.) and by Alexander von
Humboldt Foundation (I.N.), by the Spanish Consolider Ingenio 2010 Programme
CPAN (CSD2007-00042) and by MEC (Spain) under Grants FPA2007-60323,
FPA2011-23778 (O.S. and P.R.) and FPA2011-25948 (P.R.) and in part by the
funds of Polish National Science Centre under decision DEC-2011/03/B/ST2/00107
(T.P., Z.W.).

## References

  * (1) B. Aubert et al., Phys.Rev.Lett. 100 (2008) 011801; J.P. Lees et al., arXiv: 1109.1527 [hep-ex].
  * (2) M. Fujikawa et al., Phys. Rev. D78 (2008) 072006; H. Hayashii, M. Fujikawa, Nucl.Phys.Proc.Suppl. 198 (2010) 157.
  * (3) E. Braaten and C.-S. Li, Phys.Rev. D42 (1990) 3888.
  * (4) G. Aad et al., arXiv: 1207.7214 [hep-ex]; S. Chatrchyan, arXiv:1207.7235 [hep-ex].
  * (5) S. Jadach, Z. Wa̧s, R. Decker and J. H. Kühn, Comput. Phys. Commun. 76 (1993) 361.
  * (6) J.H. Kühn, E. Mirkes, Z. Phys. C56 (1992) 661; J.H. Kühn, A. Santamaría, Z. Phys. C48 (1990) 445.
  * (7) D. G. Dumm, A. Pich and J. Portolés, Phys. Rev. D69 (2004) 073002.
  * (8) T.E. Coan et al., Phys.Rev.Lett. 92 (2004) 232001.
  * (9) G. Ecker, J. Gasser, H. Leutwyler, A. Pich, E. de Rafael, Phys. Lett. B223 (1989) 425; Nucl. Phys. B321 (1989) 311.
  * (10) P. Roig, I. Nugent, T. Przedzinski, O. Shekhovtsova, Z. Was, in these Proceedings.
  * (11) O. Shekhovtsova, T. Przedzinski, P. Roig and Z. Was, arXiv: 1203.3955 [hep-ph].
  * (12) T. Przedzinski, O. Shekhovtsova and Z. Was, http://annapurna.ifj.edu.pl/$\sim$wasm/RChL/RChL.htm.
  * (13) D.G. Dumm, P. Roig, A. Pich, J. Portolés, Phys. Lett. B685 (2010) 158; Phys. Rev. D81 (2010) 034031.
  * (14) K. Nakamura, J. Phys. G37 (2010) 075021.
  * (15) I.M. Nugent, SLAC-R-936.
  * (16) E.I. Shibata et al., arXiv: hep-ex/0210039.