# One class of integrals evaluation in magnet soliton theory II.

D. Kovalenko<EMAIL_ADDRESS>University of British Columbia,
Academic 2366 Main Mall, Vancouver, BC, Canada  Zhmudsky A.A
<EMAIL_ADDRESS>Department of Physics, University of Central
Florida, 4000 Central Florida Blvd. Orlando, FL, 32816

###### Abstract

The paper proposes an approximate expression for calculating very complex one-
dimensional integrals depending on the parameter $a$. These integrals often
occur in computational problems theory of magnetic solitons.

The resulting analytical expressions are verified by numerical calculation.
The comparison shows that the relative accuracy of the approximating
expression tends to zero for large and small values of the parameter $a$. In a
small region near $a=0.119$, the error does not exceed five percent, and for
parameter values near $18.0$ the relative error does not exceed seven percent.
It turned out that the proposed method of constructing approximating
expressions can be generalized to much more complex integrals.

A program that compares an approximate expression with an exact numerical
value written in C and tested on examples that allow analytical solutions.

###### pacs:

12., 12.20.-m, 12.15.-y, 13.66.-a

##  I Introduction

Nonlinear excitations (topological solitons KIK ) play an important role in
the physics of low-dimensional magnets BI ; IvK . They contribute greatly to a
heat capacity, susceptibility, scattering cross-section and other physical
characteristics of magnets. In particular, for two-dimensional (2D) magnets
with discrete degeneracy it is important to take into account the localized
stable (with quite-long life time) 2D solitons BI ; IvK . According to the
experiments W1983 ; W1986 , these solitons determine the relaxation of
magnetic disturbance and can produce peaks in the response functions.

The traditional model describes the state of the magnet in units magnetization
vector $\vec{m}$, $\vec{m}^{2}=1$ with energy function in form

$W=\int\limits d^{2}x\left\\{A(\nabla\vec{m})^{2}+w_{0}(\vec{m})\right\\},$
(1)

where A is the inhomogeneous exchange constant, and $w_{0}(\vec{m})$ is
anisotropy energy.

In the anisotropic case, the solution is multidimensional, so the soliton
structure is determined by a system of equations in partial derivatives. There
are no general methods for finding localized solutions of such equations and
stability analysis of the solution. For this reason, the direct variational
method is often used. Consequently, the choice of trial function plays a key
role in such an analysis ZyI ; IZ ; Steph . In most cases, the choice of trial
function in the form:

$\displaystyle tg{\theta\over 2}={R\over
r}\exp\left(-\frac{r}{b}\right)(1+C_{1}\cos 2\chi),$ (2)
$\displaystyle\varphi=\chi+C_{2}\sin 2\chi+\varphi_{0},$ (3)

is quite successful. In this formula: $r$, $\chi$ are polar coordinates in the
magnetic plane, $R$, $a$, $C_{1}$, $C_{2}$ and $\varphi_{0}$ are variable
parameters.

In the papers ZyI ; IZ ; Steph  Newton’s iterative method for solving a system
of non-linear algebraic equations teuk  was used to find variable parameters
that provide a minimum of energy. The numerical algorithm indicated above
leads to the necessity of multiple calculation of two-dimensional integrals
(over the angle $\chi$ and the radius $r$). Experience shows that the
calculation process can be significantly accelerated (and the accuracy of
calculations improved) if an analytic-numerical method for calculating
integrals over the radius is prepared in advance. The typical shape of these
integrals:

$\int\limits_{0}^{\infty}{f(r)dr\over\left[(r\exp(r))^{2}+b^{2}\right]^{n}}$
(4)

where $n$ is an integer, $f(r)=\exp(r),r\exp(r),\\\ (1+r)\exp(r),r^{2}\exp(r)$
and etc. Note that $a$ is an essential parameter - there is no replacement
that could eliminate it. Below we consider the case $n=1$. On the one hand, it
is almost always possible to reduce the calculation for $n>1$ to the sum of
integrals for $n=1$. On the other side it is not difficult to modify the
calculation scheme for the case $n>1$.

In a previous work One  it was shown that with an appropriate choice of the
integration contour in the complex plane, followed by numerical calculation of
the sums of residues of the integrand, one can find values of the integral of
the type (4) with a given accuracy. However, in this case it is impossible to
write an analytic expression for the roots in the complex plane. This leads to
the need to take into account the number of poles of the integrand that is not
specified in advance.

In this paper, we propose a simple analytical expression that does not require
the calculation of a predetermined number of poles and allows to find the
value of the integral with an accuracy of no more than a few percents.

As an example, consider two integrals:

$\displaystyle
I(a)=\int\limits_{0}^{\infty}\frac{x\exp(x)dx}{(x\exp(x))^{2}+a^{2}},$ (5)
$\displaystyle
J(a)=\int\limits_{0}^{\infty}\frac{\exp(x)dx}{(x\exp(x))^{2}+a^{2}}$ (6)

It is easy to show that their sum is

$\displaystyle I(a)+J(a)=\frac{\pi}{2a}\quad$ (7)

This equality will be useful for controlling the accuracy of numerical
integration.

##  II For small values of the parameter $a$

For parameter values $a\to 0$, one can obtain a good approximation for
integrals (5,6). Let us introduce a new variable $u=x\exp(x)$. The solution of
this equation with respect to $x$ is expressed in terms of the Lambert
function $x=W_{0}(u)$ Lam ; Eul ; Dub . Integrals (5,6) become:

$\displaystyle
I(a)=\int\limits_{0}^{\infty}\frac{W_{0}(u)}{1+W_{0}(u)}\frac{du}{u^{2}+a^{2}},J(a)=\int\limits_{0}^{\infty}\frac{1}{1+W_{0}(u)}\frac{du}{u^{2}+a^{2}}$
(8)

where $W_{0}(u)$ is one of the principal branches of the Lambert function.

Let us rewrite the integrals in a form:

$\displaystyle
I(a)=\frac{\pi}{2a}\int\limits_{0}^{\infty}\frac{W_{0}(u)}{1+W_{0}(u)}\frac{2}{\pi}\frac{adu}{u^{2}+a^{2}},$
(9) $\displaystyle
J(a)=\frac{\pi}{2a}\int\limits_{0}^{\infty}\frac{1}{1+W_{0}(u)}\frac{2}{\pi}\frac{adu}{u^{2}+a^{2}}$
(10)

The following multiplier:

$\displaystyle\frac{2}{\pi}\frac{a}{u^{2}+a^{2}}\qquad$ (11)

when $a\to 0$ can be treated as delta function normalized from $0$ to
$\infty$.

So, in the limit $a\to 0$ we get:

$\displaystyle
I(a)=\frac{\pi}{2a}\left[\frac{W_{0}(a)}{1+W_{0}(a)}\right]_{a=0}=0,$ (12)
$\displaystyle
J(a)=\frac{\pi}{2a}\left[\frac{1}{1+W_{0}(a)}\right]_{a=0}=\frac{\pi}{2}$ (13)

In formulas (12,13) we take into account that $\lim\limits_{u\to
0}W_{0}(u)=u$.

Further, we show that these expressions describe well the integrals (5, 6) not
only for small values of $a$, but also for any $a>0$. Also, comparison with
the results of numerical calculations shows that the relative error of these
formulas does not exceed $7\%$ for any values of $a>0$.

##  III Approximate analytic expressions

To find approximate analytic expressions, we integrate the first of the
integrals (5) with respect to the parameter a:

$\displaystyle\int
I(a)da=\int\int\limits_{0}^{\infty}\frac{x\exp(x)dx}{(x\exp(x))^{2}+a^{2}}da=$
$\displaystyle\int\limits_{0}^{\infty}\arctg\displaystyle\left(\frac{a}{x\exp(x)}\right)dx$

The following figure shows the graph of the function
$\arctg\displaystyle\left(\frac{a}{x\exp(x)}\right)$ for parameter
$a=10,10^{2},10^{3},10^{4},10^{5}$

Function $\arctg\displaystyle\left(\frac{a}{x\exp(x)}\right)$ has a plateau
equal to $\pi/2$ for small values of the parameter $a$. Near the value
$\pi/4$, the function has an inflection point and tends to zero as
$x\to\infty$. Near the point where $x\exp(x)=a$, its argument is equal to one.

The corresponding solution for $x$ is found from the equation $x=W_{0}(a)$,
where $W_{0}(a)$ is one of the real branches of the Lambert function
$W_{0}(z)$ Lam ; Eul ; Dub .

The behavior of the integrand
$\arctg\displaystyle\left(\frac{a}{x\exp(x)}\right)$ is very close to behavior
of the Heaviside function. We replace the integrand
$\arctg\displaystyle\left(\frac{a}{x\exp(x)}\right)$ with the Heaviside
function, which is equal to $\pi/2$ for $x<W_{0}(a)$, has a jump for
$x=W_{0}(a)$ and vanishes for x greater than $W_{0}(a)$. Thus we get:

$\displaystyle\int\limits_{0}^{\infty}\arctg\displaystyle\left(\frac{a}{x\exp(x)}\right)dx\approx\frac{\pi}{2}\int\limits_{0}^{W_{0}(a)}dx=\frac{\pi}{2}{W_{0}(a)}$
(14)

After differentiation, we get:

$\displaystyle
I(a)=\frac{\pi}{2}\frac{W_{0}(a)}{a(1+W_{0}(a))}=\frac{\pi}{2a}\frac{W_{0}(a)}{1+W_{0}(a)}$
(15)

Respectively

$\displaystyle J(a)=\frac{\pi}{2a}\frac{1}{1+W_{0}(a)}$ (16)

The same result can be obtained if we calculate the area of the trapezoid
bounded by the abscissa axis, the y-axis, the horizontal line $\pi/2$ and the
tangent to the integrand at the inflection point.

In the limit $x\to 0$ expressions (15, 16) coincide with (12, 13).

##  IV Numerical calculation of integrals.

Gaussian quadrature formulas are used to calculate integrals (5,6). The limits
of integration are converted to the interval $[-1,+1]$ using the substitution:

$\displaystyle\int\limits_{0}^{+\infty}f(x)\,dx=\begin{vmatrix}x=\displaystyle\frac{1+t}{1-t}\\\
dx=\displaystyle\frac{2dt}{(1-t)^{2}}\end{vmatrix}=2\int\limits_{-1}^{+1}\frac{f\left(\displaystyle\frac{1+t}{1-t}\right)}{(1-t)^{2}}\,dt=$
$\displaystyle
2\frac{b-a}{2}\sum_{i=1}^{n}w_{i}f\left(\frac{1+\xi_{i}}{1-\xi_{i}}\right)\frac{1}{(1-\xi_{i})^{2}}$

that is, the function value must be calculated from the argument
$x=\displaystyle\frac{1+t}{1-t}$, and the resulting value must be divided by
$(1-t)^{2}$. The value of $t$ varies from $-1$ to $+1$.

The nodes and weights of the quadrature formulas are set using the
Gauss_Legendre() function (see teuk [p. 152]), which allows you to install
almost any number of integration nodes.

##  V Comparison with numerical calculations

The comparison results are presented in three tables: for small values of the
parameter $a$, values of the parameter $1\leq a\leq 10$ and for large values
of $a>10$.

In each table the first column contains the value of the parameter $a$, then
the numerical value of the integral, then the approximation by the formula
(16) and the relative error.

The relative error in this table is calculated as the difference between the
approximate expression and the numerical value, referring to the numerical
value. That is, the minus sign indicates that the numerical value is greater
than the approximating expression.

###  V.1 For small values of the parameter $a$

For small values of the parameter $a$ ($0<a<1.0$), quite acceptable results
were obtained:

parameter | Exact | Approximate | Relative
---|---|---|---
a | value | value | accuracy (%)
$1\cdot 10^{-9}$ | 1.5707963066 | 1.5707963252 | 1.18e-06
$1\cdot 10^{-8}$ | 1.5707961484 | 1.5707963111 | 1.03e-05
$1\cdot 10^{-7}$ | 1.5707947727 | 1.5707961697 | 8.89e-05
$1\cdot 10^{-6}$ | 1.5707830885 | 1.5707947560 | 7.42e-04
$1\cdot 10^{-5}$ | 1.5706869694 | 1.5707806191 | 5.96e-03
$1\cdot 10^{-4}$ | 1.5699329829 | 1.5706392786 | 4.49e-02
$5\cdot 10^{-3}$ | 1.5644626729 | 1.5692286650 | 3.04e-01
$1\cdot 10^{-2}$ | 1.5302197122 | 1.5553956178 | 1.64 %
$0.119$ | 1.3746503969 | 1.4394118152 | 4.71 %
$1.0$ | 1.0031969141 | 1.0023310162 | -8.63e-02 %

The relative error tends to zero as $a\to 0$.

The largest error value does not exceed five percent (at
$a\leavevmode\nobreak\ =\leavevmode\nobreak\ 0.119$).

###  V.2 For the parameter value $1\leq a\leq 10$

For the values of parameter $a$ within ($1\leq a\leq 10$), the results
obtained are presented in the table below.

parameter | Exact | Approximate | Relative
---|---|---|---
a | value | value | accuracy (%)
1.00 | 1.0031969141e+00 | 1.0023310162e+00 | -8.63e-02
2.00 | 8.7427271092e-01 | 8.4788495181e-01 | -3.01e+00
3.00 | 8.0223081578e-01 | 7.6627616508e-01 | -4.48e+00
4.00 | 7.5357260455e-01 | 7.1329545123e-01 | -5.34e+00
5.00 | 7.1747734826e-01 | 6.7511053210e-01 | -5.90e+00
6.00 | 6.8913506420e-01 | 6.4577916569e-01 | -6.29e+00
7.00 | 6.6600973973e-01 | 6.2225892231e-01 | -6.56e+00
8.00 | 6.4661048691e-01 | 6.0280493336e-01 | -6.77e+00 %
9.00 | 6.2999141037e-01 | 5.8633322110e-01 | -6.92e+00 %
10.0 | 6.1551743151e-01 | 5.7212904958e-01 | -7.04e+00

The largest absolute value of the relative error is seven percent.

All marked results were obtained using the formula (16).

###  V.3 For large values of $a$

The table shows the values of the relative error (in percent) for the
parameter values $a$ exceeding $10^{2}$.

parameter | Exact | Approximate | Relative
---|---|---|---
a | value | value | accuracy (%)
$1\cdot 10^{1}$ | 6.1551743151e-01 | 5.7212904958e-01 | -7.04
$1\cdot 10^{2}$ | 3.83403357e-01 | 3.5816890083e-01 | -6.58
$1\cdot 10^{3}$ | 2.62958980e-01 | 2.5134338355e-01 | -4.41
$1\cdot 10^{4}$ | 1.96425062e-01 | 1.9081944919e-01 | -2.85
$1\cdot 10^{5}$ | 1.55727687e-01 | 1.5273327991e-01 | -1.92
$1\cdot 10^{6}$ | 1.28603919e-01 | 1.2684736369e-01 | -1.36
$1\cdot 10^{7}$ | 1.09333555e-01 | 1.0822372170e-01 | -1.01
$1\cdot 10^{8}$ | 9.49773498e-02 | 9.4234605334e-02 | -7.82e-01
$1\cdot 10^{9}$ | 8.38880237e-02 | 8.3367963716e-02 | -6.19e-01
$1\cdot 10^{10}$ | 7.50753831e-02 | 7.4697789991e-02 | -5.02e-01
$1\cdot 10^{11}$ | 6.79100965e-02 | 6.7627675830e-02 | -4.15e-01

The relative error in this table is calculated as the difference between the
approximate expression and the numeric value referring to the numeric value.
The minus sign indicates that the numeric value is greater than the
approximating expression.

In the picture below, you can see the relative error behavior on a logarithmic
scale. Near $a=0.1$, its maximum value is less than $5\%$. And near $a\approx
22.0$ the maximum error is less than $7\%$ in absolute value. For $a\to 0$ and
$a\to\infty$ the relative error tends to zero according to (16) and (15).

##  VI Program description

The text of the program consists of the main function and eight functions:

  * •

Gauss_Legendre \- finds the nodes and weights of the Gauss quadrature formula.
In this case, the number of nodes is 24.

  * •

gaussGL \- calculates the integral sum over the given nodes and weights.

  * •

sect \- divides the integration segment to satisfy the precision of the
integral calculation.

  * •

topoI \- the function calculates the integrand (5).

  * •

topoJ \- the function calculates the integrand (6).

  * •

bisection \- finds the value of the Lambert function by the value of the
argument $x=W_{0}(a)$.

  * •

Lambert \- calculates the difference $x\exp(x)-a$. Used by the bisection
function, which looks for the root of the equation $x\exp(x)=a$.

  * •

sign is a helper function for bisection. Specifies the sign of the difference
$x\exp(x)-a$.

When testing the program, the same trial functions were used as in the
previous paper One .

###  VI.1 Main function

In the text of the main function, all auxiliary print statements are removed.
The results of the calculations are printed to the LambertTopo.txt file.

    int main()
    {  Node = (double*)calloc(NodeNumber, sizeof(double));
       Weight = (double*)calloc(NodeNumber, sizeof(double));
       Gauss_Legendre(-1.0,1.0,Node,Weight,NodeNumber);
       double Xmin = -1.0, rel = 0.0;
       double Xmax = 1.0;
       irecmax = nrec = 60;
       QuadratureRule = gaussGL;
       double int0;
       for(int m = 0; m < nPoints; m++)
        {  for(int k = 0; k < 2; k++)
           {  f1 = fun[k];
        Xmin = Left[k];
        Xmax = Right[k];
        estimation = fabs(QuadratureRule(Xmin,Xmax));
        int0 = sect(Xmin,Xmax,estimation);
        topoIJ[k][m] = A*int0;
        if(m > 0) rel = (topoIJ[k][m] - topoIJ[k][m+1])/M_LN10;
        nofun = unworked = 0;
        irec = 0;
        double  a = 0.0, c = 25.0;
        double funa, func, root;
        funa = Lambert(a);
        func = Lambert(c);
        root = W0 = bisection(a,c,funa,func,Lambert);
        W[m] = root;
        error[m] = 100*(M_PI_2/(1+W[m]) - topoIJ[1][m])/topoIJ[1][m];
        A += 0.001;
        if(m == 0) fprintf(ass," topoIJ[0]=%8.3lf  topoIJ[1]=%8.3lf  \n",
                    topoIJ[0][m],topoIJ[1][m]);
          }
        A = Astart;
        fprintf(ass,"* topoIJ[0]=%8.3lf  topoIJ[1]=%8.3le  \n",
                    topoIJ[0][0],topoIJ[1][0]);
        fprintf(ass,"       A            topoI[m]        M_PI_2 - topoI[m]
                    topoJ[m]       M_PI_2/(1+W(a))           W(a)
                    error[m] (%%)     topoI[m]+topoJ[m]  \n");
        for(int m = 0; m < nPoints; m++)
        {fprintf(ass," %8.3le   %18.8le   %18.8le  %20.10le %20.10le
                    %20.10le %20.10le %20.10le\n",A,topoIJ[0][m],
                    M_PI_2 - topoIJ[0][m],topoIJ[1][m],M_PI_2/(1+W[m]),
                    W[m],error[m],topoIJ[0][m]+topoIJ[1][m]);
            A += 0.001;
        }
        fclose(ass);
        free(Node);
        free(Weight);
        return 1;
    }

During the program execution: we set the values of the nodes and weights of
the Gaussian quadrature formulas, find the values of the integrals, calculate
the value of the Lambert function for a given $a$, calculate the approximating
expression and the error value (in percent), and print the results.

###  VI.2 Recursive adaptive quadrature algorithm

The integration algorithm in the program consists of two independent parts:
the adaptive function (sect(double, double, double)) and the quadrature
formula (gaussGL(double, double)). The adaptive function uses a recursive
algorithm to implement the standard bisection method. This algorithm has been
tested many times in previous papers One ; BZKINR ; BZ ; BZh .

To achieve the required relative accuracy $\varepsilon$ of integration, the
integral estimate $\sum_{whole}$ found on the interval [$a_{i}$,$b_{i}$] is
compared with the sum of $\sum_{left}$ and $\sum_{right}$, which are
calculated in the left and right parts of the segment [$a_{i}$,$b_{i}$]. The
comparison rule was chosen in the form:

$\sum_{left}+\sum_{right}-\sum_{whole}\leq\varepsilon\cdot\sum_{whole},$

where $\sum_{whole}$ denotes the integral sum over the entire integration
interval [$a_{i}$,$b_{i}$]. The value of $\sum_{whole}$ is accumulated and
adjusted at each halving step.

To calculate the integral sums, we use the recursive algorithm for binary
division of a segment in half, which has already been used in previous works
BZKINR ; BZ ; BZh :

    double sect(double x, double X, double int0)
    { if(++irec < irecmax)
    { double sec = (x+X)/2;
        double int1 = QuadratureRule(x,sec);
        double int2 = QuadratureRule(sec,X);
        if((fabs(int1+int2-int0) <= EPS*estimation)) int0 = int1 + int2;
        else int0 = sect(x,sec,int1) + sect(sec,X,int2);
    }
    else unworked++;
        irec--;
        return int0;
    }

The choice is due to the simplicity of this function and good test results on
numerous examples.

###  VI.3 Lambert function evaluation

To calculate the values of the Lambert function:

$x\exp(x)=u,\qquad\Rightarrow\qquad x=W_{0}(u)$

the bisection algorithm, designed as a recursive function, was used

    double bisection(double& left, double& rigth,
             double& fun_left, double& fun_rigth, double (*fun)(double&))
    { center = 0.5*(left + rigth);
        if(++irec < irecmax)
        {if(fabs(left - rigth) < 1.0e-14*fabs(center)) return center;
            double fun_center = fun(center);
            if(fabs(fun_left - fun_rigth) < 1.0e-14*fabs(fun_center))
                    return center;
            if(sign(fun_left) == sign(fun_center))
                    { left = center ; fun_left = fun_center;  }
            else { rigth = center;  fun_rigth = fun_center;  }
            bisection(left,rigth,fun_left,fun_rigth,fun);
        }
        irec--;
        return center;
    }
    double Lambert1(double& x)
    {     return x*exp(x) - A;   }

##  VII More generali integrals

It is easy to see that calculations carried out in Sections II and III can be
extended to a wider class of integrals. Consider the following example:

$\displaystyle
K(a)=\frac{1}{a}\int\limits_{0}^{\infty}g(x)\frac{a}{f^{2}(x)+a^{2}}dx$ (17)

where $f(x)$ is a continuous function that increases monotonically from zero
to $\infty$ as the value of $x$ changes within the same limits. We also assume
that function $g(x)$ does not break the convergence of the integral. Let us
introduce the new variable $u=f(x)$ and denote the solution of the equation
$f(x)=u$ as $F(u)$, where $F(u)$ is the inverse function of $f(x)$.

After substitution into (17) we get:

$\displaystyle
K(a)=\frac{1}{a}\int\limits_{0}^{\infty}g(F(u))\frac{a}{u^{2}+a^{2}}\frac{du}{f^{\prime}(F(u))}$
(18)

Let’s rewrite the integral as follows:

$\displaystyle
K(a)=\frac{\pi}{2a}\int\limits_{0}^{\infty}\frac{g(F(u))}{f^{\prime}(F(u))}\frac{2}{\pi}\frac{adu}{u^{2}+a^{2}}$
(19)

###  VII.1 Asymptotic value at $a\to 0$

When $a\to 0$ the following multiplier:

$\displaystyle\frac{2}{\pi}\frac{a}{u^{2}+a^{2}}\qquad$ (20)

can be considered as delta function normalized from $0$ to $\infty$. Which
results in:

$\displaystyle
K(a)=\frac{\pi}{2a}\left[\frac{g(F(u))}{f^{\prime}(F(u))}\right]_{a=0}$ (21)

In the particular case of $f(x)=x\exp(x)$ and $g(x)=x\exp(x)$ from this
expression yields the result obtained in section II.

###  VII.2 Approximate expression

Under unchanged conditions on the functions $f(x)$ and $g(x)$, consider the
same integral:

$\displaystyle K(a)=\int\limits_{0}^{\infty}\frac{g(x)}{f^{2}(x)+a^{2}}dx$
(22)

Following the same calculations as in the III section, we integrate (22) with
respect to the parameter a:

$\displaystyle\int
K(a)da=\int\int\limits_{0}^{\infty}\frac{g(x)}{f(x)}\frac{dx}{1+\left[\displaystyle\frac{a}{f(x)}\right]^{2}}d\left[\frac{a}{g(x)}\right]=$
$\displaystyle\int\limits_{0}^{\infty}\left[\frac{g(x)}{f(x)}\right]\arctg\displaystyle\left(\frac{a}{f(x)}\right)dx$

Since the function $f(x)$ increases monotonically from zero to $\infty$ as the
value of $x$ changes within the same limits, the behavior of function
$\arctg\displaystyle\frac{a}{f(x)}$ is the same as in the picture in the
section III.

Therefore, we can approximately replace $\arctg\displaystyle\frac{a}{f(x)}$
with the Heaviside function. The iflection point satisfies the equation
$a=f(x)$.

This equation can be solved with respect to $x$ as $x=F(a)$, where $F(a)$ is
the inverse function to the $f(x)$. So we can integrate from 0 to $F(a)$, and
the integral reduces to:

$\displaystyle\int
K(a)da=\frac{\pi}{2}\int\limits_{0}^{F(a)}\left[\frac{g(x)}{f(x)}\right]dx$

where $F(a)$ is the inverse function of $a=f(x)$. Using the integral
differentiation formula with the upper limit depending on the parameter, we
find:

$\displaystyle
K(a)=\frac{\pi}{2}\left[\frac{g(x)}{f(x)}\right]_{F(a)}\cdot\frac{dF(a)}{da}$

Let’s use the obtained formula to calculate the integral (14). In this case
$F(a)=W_{0}(a)$ and $g(x)=x\exp(x)$. The derivative of the inverse function is
equal to one divided by the derivative of $da/dx$ after $x$ has been replaced
by $a$. Thus:

$\displaystyle\frac{da}{dx}=(1+x)\exp(x)\quad\Rightarrow\quad\frac{dx}{da}=\frac{1}{(1+a)\exp(a)}$
$\displaystyle K(a)=\frac{\pi}{2a}\frac{W_{0}(a)}{a(1+W_{0}(a))}$

Which is the same as the previously obtained expression (15).

##  VIII Conclusion

The method of finding the approximating expression proposed in the sections II
and III can be used for not only for (5, 6) integrals, but also for many other
similar integrals. For example, all integrals of this type that have a
denominator higher than the first order:

$\displaystyle
I(a)=\int\limits_{0}^{\infty}\frac{x\exp(x)dx}{[(x\exp(x))^{2}+a^{2}]^{n}},$
(23) $\displaystyle
J(a)=\int\limits_{0}^{\infty}\frac{\exp(x)dx}{[(x\exp(x))^{2}+a^{2}]^{m}}$
(24)

can be represented as derivatives of the appropriate order with respect to the
parameter $a$.

Also, the first term in the denominator does not necessarily have to be
$x\exp(x)$. Obviously, this can be any continuous function $f(x)$
monotonically increasing from zero to $\infty$ while the value of $x$ changes
within the same limits. For example, $\sinh(x)$, $x\cosh(x)$ etc.

##  IX Acknowledgments

We expresse our thanks to Dr. V.K.Basenko for stimulating and usefull
discussion of the problem.

## References

  * (1) A.M.Kosevich, B.A.Ivanov, and A.S.Kovalev.  Magnetic solitons.  Phys. Rep, 194(3-4):117 – 238, 1990.
  * (2) V.G. Baryakhtar and B.A. Ivanov.  Soliton thermodynamics of low-dimensional magnets.  Soviet Scintific Review.Section A Phys. Rev.Ed. by I.Khalatnikov., 16(3):192, 1992.
  * (3) B.A. Ivanov and A.M. Kosevich.  Solitons in low-dimensional antiferromagnets.  Fizika Nizkikh Temperatur, 21:355 – 389, 1995.
  * (4) F Waldner.  Two-dimensional soliton contribution to the esr linewidth in layered antiferromagnets.  Journal of Magnetism and Magnetic Materials, 31:1203 – 1204, 1983\.
  * (5) F Waldner.  Two dimensional soliton energy and esr in afm.  Journal of Magnetism and Magnetic Materials, 54:873–874, 1986.
  * (6) A. A. Zhmudskii and Boris A. Ivanov.  On the structure and stability of two-dimensional dynamic solitons in ferromagnets.  Journal of Experimental and Theoretical Physics Letters, 65:945–950, 1986.
  * (7) A. Zhmudskii and B. Ivanov.  On the structure and stability of two-dimensional dynamic solitons in ferromagnets.  JETP Letters, 65(6):945–950, 1986.
  * (8) B.A. Ivanov, V.A. Stephanovich, and A.A. Zhmudskii.  Magnetic vortices. the microscopic analogs of magnetic bubbles.  Journal of Magnetism and Magnetic Materials, 88(1):116–120, 1990\.
  * (9) W.H. Press, S.A. Teukolsky, W.T. Vetterling, and B.P. Flannery.  Numerical Recipes in C.  Cambridge University Press, Cambridge, 1992.
  * (10) A.A. Zhmudskii.  One class of integrals evaluation in magnet solitons theory.  https://arxiv.org/abs/physics/981, (9812003v1):122, December 2 1998\.
  * (11) Lambert J. H.  Observationes variae in mathesin puram.  Acta Helveticae physico-mathematico-anatomico-botanico-medica, 3:128–168, 1758.
  * (12) L. Euler.  De serie lambertina plurimisque eius insignibus proprietatibus.  Acta Academiae Scientarum Imperialis Petropolitinae, 2:2951, 1773\.
  * (13) A. E. Dubinov, I. D. Dubinova, and Saykov S. K.  The Lambert W Function and Its Applications to Mathematical Problems of Physics (in Russian.  RFNC-VNIIEF, Sarov, Russia, 2006.
  * (14) A.N Berlizov and A.A.Zhmudsky.  The recurcive one-dimensional adaptive quadrature code.  Preprint KINR-98-4, Kiev Institute for Nuclear Reasearch,, 1998.
  * (15) A.N Berlizov and A.A.Zhmudsky.  The recursive adaptive quadratures in fortran-77.  Athens, GREECE, July 4-8 1999. ggg.  3 rd International Conference on COMPUTERS.
  * (16) A.N Berlizov and A.A.Zhmudsky.  The recursive adaptive quadrature in ms fortran-77.  Preprint LANL, (physics/arXiv:physics/9905035), 19 May 1999.
  * (17) D.D. MacCracken and W.S Dorn.  Numerical methods and Fortran programming : with applications in engineering and science.  Wiley, New York, 1965.


aainstitutetext: Faculty of Physics, University of Warsaw, ul. Pasteura 5,
02-093 Warsaw, Polandbbinstitutetext: The Laboratory for Quantum Gravity &
Strings, Department of Mathematics & Applied Mathematics, University of Cape
Town, Cape Town, South Africaccinstitutetext: The National Institute for
Theoretical and Computational Sciences, Private Bag X1, Matieland, South
Africa

# Spread Complexity and Topological Transitions in the Kitaev Chain

Pawel Caputa b Nitin Gupta b,c S. Shajidul Haque a Sinong Liu b,c Jeff Murugan
b,c and Hendrik J.R. Van Zyl

###### Abstract

A number of recent works have argued that quantum complexity, a well-known
concept in computer science that has re-emerged recently in the context of the
physics of black holes, may be used as an efficient probe of novel phenomena
such as quantum chaos and even quantum phase transitions. In this article, we
provide further support for the latter, using a 1-dimensional p-wave
superconductor - the Kitaev chain - as a prototype of a system displaying a
topological phase transition. The Hamiltonian of the Kitaev chain manifests
two gapped phases of matter with fermion parity symmetry; a trivial strongly-
coupled phase and a topologically non-trivial, weakly-coupled phase with
Majorana zero-modes. We show that Krylov-complexity (or, more precisely, the
associated spread-complexity) is able to distinguish between the two and
provides a diagnostic of the quantum critical point that separates them. We
also comment on some possible ambiguity in the existing literature on the
sensitivity of different measures of complexity to topological phase
transitions.

##  1 Introduction and motivation

Nature is complex. Not complicated, like an aeroplane; complex, like a single
human cell. This is why much of our time as theoretical physicists is spent
searching for toy models; simple enough to solve but rich enough to contain
some universal lessons about Nature. Not all toy models are created equal.
Poor ones are simply trivial. Better ones answer the questions we set out to
answer, but the best ones are those, like the harmonic oscillator or the Ising
model, whose lessons go far beyond the problems they were conceived for. The
Kitaev chain, a 1-dimensional tight-binding chain of spinless electrons
supplemented with a superconducting term, falls into this last class.
Originally introduced by Kitaev Kitaev_2001  in order to furnish free Majorana
fermions that would facilitate topological quantum computing logic gates, this
model has since matured into a versatile laboratory in which to test all
manner of questions involving topological phases, critical points and quantum
matter more generally Leijnse:2012tn .

The Kitaev Hamiltonian, $H_{\mathrm{Kitaev}}(J,\mu,\Delta)$, is parameterized
by three real parameters; a tight-binding coupling $J$ that controls the
nearest-neighbour hopping between sites on the 1-dimensional lattice, a
chemical potential, $\mu$, and a p-wave pairing gap, $\Delta$ controlling the
superconducting term. The Hamiltonian exhibits two distinct gapped phases;
when $|\mu|>|J|$ the system finds itself in a trivial phase with a unique
ground state independent of the boundary conditions, while in its topological
phase, characterised by $|\mu|<|J|$, it manifests Majorana zero-modes in the
open Kitaev chain, and odd fermion parity of the groundstate in the closed
chain with periodic boundaries. For different values of its parameters, the
Kitaev chain shares many properties with several other well-known quantum
systems:

  * •

When $J,\Delta>0$, the Kitaev chain belongs to the same universality class as
the Ising model. More precisely, the Kitaev and transverse Ising chains can be
mapped into each other through a Jordan-Wigner transform Greiter_2014 .

  * •

The Su-Schrieffer-Heeger (SSH) model of polyacetylene Su:1979ua , is a
1-dimensional dimerized lattice of electrons in which each unit cell consists
of two sites $A$ and $B$ and a hopping term connecting the two sub-lattices.
Performing a particle-hole transformation, $c_{B,n}\to c_{B,n}^{\dagger}$ on
all electrons living on the $B$ sublattice, converts it, for specific
parameter values, into the p-wave superconductor, and

  * •

Turning off the superconducting term ($\Delta=0$) $H_{\mathrm{Kitaev}}\to
H_{\mathrm{XX}}$, the isotropic limit of the XY model Franchini:2016cxs .

In most cases these, and other, properties derive from a study of the spectrum
of eigenvalues, either analytically or numerically, of the Hamiltonian matrix.
However, the explicit diagonalization of an $N\times N$ matrix is a
computationally expensive exercise for large $N$. It behoves us then to search
for more bespoke tools to study such many-body systems, especially in the
presence of further complications such as interactions, disorder and large
system sizes. One particularly exciting set of such tools comes from the
growing influence of quantum information theory in both condensed matter and
high energy physics. Included in this set are concepts such as Von Neumann
entropy, entanglement spectra, Holevo entropy, quantum error correction and,
of particular interest to us, quantum complexity nielsen00 ; Watrous_2008 .

In a nutshell, quantum complexity - the quantum analog of the well-known
computer science concept of computational complexity - attempts to quantify,
as efficiently as possible, how difficult it is to synthesize a particular
quantum state using a set of allowed unitary operations. A substantial amount
of the recent interest in quantum complexity in high energy physics is
motivated by developments in the understanding of the physics of black hole
information loss. In particular, the AdS/CFT correspondence Maldacena:1997re
that connects quantum gravity in $D$-dimensional anti-de Sitter space to a
conformal field theory living on its $(D-1)$-dimensional boundary, also
connects quantum complexity in the field theory to various quantities in the
gravity side of the correspondence via the so-called complexity/volume and
complexity/action conjectures Brown:2015bva . Indeed, complexity as a resource
in quantum information Halpern:2021ufe , has featured prominently in
discussions on black hole firewall Harlow:2013tf , fast scrambling near black
hole horizons Aaronson:2016vto ; Susskind:2014rva ; Stanford:2014jda ;
Brown:2015bva , quantum chaos Maldacena:2015waa  and the properties of the
holographic dictionary itself Bouland:2019pvu .

Forgetting about the connection to black holes and gravitational physics for
the moment, one, not insignificant problem is that, while the central question
of finding the smallest number of unitary operations that connect a fixed
reference and target state remains, the actual computation of quantum
complexity is not unique. There are by now several contenders (see e.g.
reviews Chapman:2021jbh ; Chen:2021lnq  and references therein). Among these,
the Nielsen approach Nielsen1 ; Nielsen_2006 ; Nielsen3  phrases the problem
as one of finding the minimal geodesics connecting two points on the manifold
of quantum states. The inputs required for its computation are the desired
reference and target state along with the choices of unitary gates, their
penalties and a cost function. These ambiguities make it difficult to
determine definitively which features of a quantum system the complexity is
sensitive to. That said, Nielsen complexity has been shown already to be a
useful probe of topological phase transitions in Ali:2018aon ; Liu:2019aji ;
Xiong:2019xoh  with some variation between the different studies, possibly due
to the ambiguities outlined above111It must be pointed out though that the
authors of Liu:2019aji ; Xiong:2019xoh  have argued that appropriate choices
for the ambiguities may detect topological phase transitions as
discontinuities in Nielsen complexity.

More recently, a proper subset of us (SL & PC) have argued in Caputa:2022eye
that a recently introduced measure of state complexity called spread
complexity Balasubramanian:2022tpr  is able to disambiguate between the
topologically distinct phases of the SSH model Su:1979ua . Unlike Nielsen
complexity however, spread complexity requires as input only the desired
reference and target state along with the Hamiltonian of the system in
question making it less susceptible to the ambiguities of Nielsen complexity.
As a variant of Krylov complexity, it also comes equipped with a clear
algorithmic construction, making it computationally efficient. Taken together,
these properties make it a useful computational tool to study various aspects
of many-body quantum systems. For example, in addition to successfully probing
the topological phase structure of the SSH model, spread complexity has also
been been shown to be a useful diagnostic of quantum chaos
Balasubramanian:2022tpr  to supplement other diagnostics Ali:2019zcj ;
Bhattacharyya:2019txx ; Bhattacharyya:2020art ; Bhattacharyya:2020iic  such as
the out-of-time-order correlators (OTOCs) Larkin1969QuasiclassicalMI  or the
spectral form factor (SFF) Dyer:2016pou ; Liu:2018hlr . This article is an
extension of these ideas and provides further substantiation of spread
complexity as a good probe of quantum phase transitions.

The rest of this article is organised as follows: Since the spread complexity
that is the focus of this paper is closely related to K(rylov)-complexity, we
begin the next section with a brief self-contained review of spread complexity
based on the Krylov basis Balasubramanian:2022tpr . Having reviewed the
technology, in section 3 we diagonalize the Kitaev Hamiltonian in terms of
$SU(2)$ coherent states and subsequently compute the spread complexity of the
1-dimensional Kitaev chain. This is followed by a detailed study of the spread
complexity of the p-wave superconductor for three different circuits,
involving different choices of target state. Finally we conclude and discuss
future directions. Appendix A gives a self-contained treatment of the $SU(2)$
coherent states Perelomov , Gazeau:2009zz  utilized in the computation of the
spread complexity. Finally, in Appendix B we apply the methods developed in
Caputa:2022eye  to a generalization of the Kitaev chain discussed above, with
the open/closed boundary conditions replaced by twisted boundary conditions
Kawabata:2017zsb . The phase that enters the problem with the twisting,
together with the p-wave pairing gap parameter $\Delta$ endow the Kitaev chain
further structure beyond that of the SSH model. We provide some preliminary
evidence that some of this structure manifests in the spread complexity.

##  2 Krylov and Spread Complexity

In the interests of self-containment, we provide here a brief summary of the
notion of spread complexity introduced in Balasubramanian:2022tpr  as well as
the Krylov basis, both of which are central to the computations performed in
this paper. This section borrows heavily from a similar discussion in
Caputa:2022eye .

Our starting point is a general quantum state $|\Psi(s)\rangle$, related to
some initial state $|\Psi_{0}\rangle$ by a unitary transformation

$|\Psi(s)\rangle=e^{-isH}|\Psi_{0}\rangle\,.$ (1)

We will refer to the state $|\Psi_{0}\rangle$ as the reference state and
$|\Psi(s)\rangle$ as the target state of a unitary circuit with circuit
Hamiltonian, $H$. The parameter $s$ is a circuit parameter that we will
usually take to run from $0$ to $1$ as we traverse the circuit. A useful
measure of quantum state complexity can be defined as the spread of
$|\Psi(s)\rangle$ in the Hilbert space Balasubramanian:2022tpr , i.e., the
spread of $|\Psi_{0}\rangle$ under the unitary evolution (1). To quantify this
spread, the target state may be expanded in some basis
$\mathcal{B}=\left\\{|B_{n}\rangle,n=0,1,2,\cdots\right\\}$ where we take
$|B_{0}\rangle=|\Psi_{0}\rangle$, the reference state. Then, we can define a
cost function,

$C_{\mathcal{B}}(s)=\sum_{n}n|\langle\Psi(s)|B_{n}\rangle|^{2}\,,$ (2)

which is to be minimized over all choices of basis $\mathcal{B}$. We then
define the spread complexity as

$\mathcal{C}(s)=\textnormal{min}_{\mathcal{B}}\ C_{\mathcal{B}}(s),$ (3)

and the actual value of complexity corresponding to our target state is
$\mathcal{C}(s=1)$.

The basis that satisfies the minimisation condition of the cost function (2)
is the Krylov basis Balasubramanian:2022tpr . The Krylov basis
$\\{|K_{n}\rangle\\}$ can in turn be obtained by performing a Gram-Schmidt
orthogonalization on the states $|O_{n}\rangle=H^{n}|\Psi_{0}\rangle$. This is
known as the Lanczos algorithm Lanczos1950AnIM  (see review
viswanath2008recursion ). Note that the Krylov basis depends on the choice of
circuit Hamiltonian as well as the choice of reference state. The circuit
Hamiltonian is generally tri-diagonal in this basis in the sense that

$H|K_{n}\rangle=a_{n}|K_{n}\rangle+b_{n}|K_{n-1}\rangle+b_{n+1}|K_{n+1}\rangle\,,$
(4)

where the so-called Lanczos coefficients $a_{n}$ and $b_{n}$ may be
recursively constructed using the Lanczos algorithm. The information captured
by these coefficients is also contained in the moments of the return-amplitude

$S(s)\equiv\langle\Psi(s)|\Psi(0)\rangle=\langle\Psi_{0}|e^{isH}|\Psi_{0}\rangle,$
(5)

that is the probability amplitude for the state to remain in the initial
state. In many-body literature it is also often called the survival or
Loschmidt amplitude Heyl:2017blm .

Then, the target state may be expanded in the Krylov basis

$|\Psi(s)\rangle=\sum_{n}\psi_{n}(s)|K_{n}\rangle\,,$ (6)

where, by construction, the coefficients satisfy the discrete Schrodinger
equation

$i\partial_{s}\psi_{n}(s)=a_{n}\psi_{n}(s)+b_{n}\psi_{n-1}(s)+b_{n+1}\psi_{n+1}(s)\,.$
(7)

With the knowledge of the Lanczos coefficients this equation may be solved
with initial condition $\psi_{n}(0)=\delta_{n,0}$. Finally, in the Krylov
basis the spread complexity now becomes

$\displaystyle\mathcal{C}(s)$ $\displaystyle=$
$\displaystyle\sum_{n}n|\psi_{n}(s)|^{2}\,.$ (8)

This definition generalises Krylov complexity (K-complexity for short) of the
operator growth defined in Parker:2018yvk  (see also Barbon:2019wsy ;
Magan:2020iac ; Rabinovici:2020ryf ; Jian:2020qpp ; Yin:2020oze ;
Dymarsky:2019elm ; Dymarsky:2021bjq ; Caputa:2021ori ; Carrega:2020jrk ;
Kim:2021okd ; Kar:2021nbm ; Rabinovici:2022beu ; Trigueros:2021rwj ;
Rabinovici:2021qqt ; Hornedal:2022pkc ; Bhattacharjee:2022vlt ; Fan:2022xaa ;
Adhikari:2022whf ; Muck:2022xfc ; Patramanis:2021lkx ; Caputa:2021sib ;
Bhattacharjee:2022qjw  and references therein) 222There, by appropriate choice
of the inner-product in the space of operators one can set $a_{n}=0$.. It is
often useful to introduce a "complexity operator" by relations

$\hat{K}=\sum_{n}n|K_{n}\rangle\langle
K_{n}|,\qquad\mathcal{C}(s)\equiv\langle\Psi(s)|\hat{K}|\Psi(s)\rangle\,.$ (9)

These two relations ((8) and (9)) will be our working definitions of the
spread complexity for the remaining parts of this work.

In principle, every quantum dynamics can be formulated in terms of this one-
dimensional chain-like picture. In practice, for generic setups, extracting
Lanczos coefficients and solving the Schrodinger equation may be only
accessible numerically. Nevertheless, as it will also be the case in this
work, when the circuit Hamiltonian is an element of a low-rank semi-simple Lie
algebra, the Krylov basis and complexity may be straightforwardly extracted
following the construction reviewed in Caputa:2021sib . In Appendix (A) we
also show how the spread complexity may be computed directly from the relevant
coherent state overlap.

Last but not the least, as in Caputa:2022eye , we will employ a natural
extension of the above definitions to free theories. Namely, the states of our
interest will have the product form for each momentum mode. Similarly, to
geometric definitions (Nielsen type) we will then compute the spread
complexity for each momentum separately and then just sum the answers (or
integrate in the continuum limit)333We should admit that, for general states,
demonstrating that this prescription is equivalent to the canonical definition
of K-complexity in position space or the definition starting from the return
amplitude is not immediately obvious. In this work we will then follow the
steps outlined above and leave studies of relations between different
prescriptions as an important and interesting future direction. .

##  3 The Kitaev chain and topological phase transitions

In this section we review our main physical setup of the Kitaev model
Kitaev_2001 . The Kitaev chain Hamiltonian is given by

$H=\sum_{j=1}^{L}\left[-\frac{J}{2}\left(c_{j}^{\dagger}c_{j+1}+c_{j+1}^{\dagger}c_{j}\right)-\mu\left(c_{j}^{\dagger}c_{j}-\frac{1}{2}\right)+\frac{1}{2}\left(\Delta
c_{j}^{\dagger}c_{j+1}^{\dagger}+\Delta^{*}c_{j+1}c_{j}\right)\right]\,,$ (10)

where the $c_{i}$ ($c_{i}^{\dagger}$) are spinless fermionic annihilation
(creation) operators acting at site $i$ on an $L$-site lattice, defined by the
algebra,

$\left\\{c_{j}^{\dagger},c_{k}\right\\}=\delta_{jk}.$ (11)

The hopping amplitude, $J$ and chemical potential $\mu$ are real parameters
while the superconducting pairing strength, $\Delta$, is generally complex,
although it can always be chosen to be a positive real number by rescaling the
annihilation operator $c_{i}\to
e^{i\frac{1}{2}\operatorname{arg}\Delta}c_{i}$. In what follows, we consider
$\Delta$ real which is equivalent to a rescaling factor of either
$e^{i\frac{1}{2}\operatorname{arg}\Delta}=\pm 1$ or
$e^{i\frac{1}{2}\operatorname{arg}\Delta}=\pm i$, for positive or negative
$\Delta$, respectively. The model exhibits a topological phase transition when
the quasi-particle spectrum is gapless Kitaev_2001 . In the continuum limit,
$L\rightarrow\infty$, this corresponds to a topological phase transition at
$|\frac{\mu}{J}|=1$ with a topological phase in the regime $|\frac{\mu}{J}|<1$
and a non-topological phase in the regime $|\frac{\mu}{J}|>1$. Furthermore,
the winding number of the topological phase is determined by the sign of
$\Delta$, which implies that there are two topological phases in the
$\frac{\Delta}{J}-\frac{\mu}{J}$ plane, separated by $\Delta=0$.

To facilitate a discussion of the spectral properties of the model, it is
convenient to perform a Fourier transform

$c_{j}=\frac{1}{\sqrt{L}}\sum_{n}e^{ik_{n}j}a_{k_{n}}\,,$ (12)

which puts the Hamiltonian into the form

$\begin{split}H=&\sum_{n=0}^{L-1}\left[-(\mu+J\cos(k_{n}))\left(a^{\dagger}_{k_{n}}a_{k_{n}}-\frac{1}{2}\right)+\frac{1}{2}\Delta\left(e^{-ik_{n}}a^{\dagger}_{-k_{n}}a^{\dagger}_{k_{n}}+e^{ik_{n}}a_{k_{n}}a_{-k_{n}}\right)\right]\\\
=&\sum_{n=0}^{\lceil\frac{L-1}{2}\rceil}\left[-(\mu+J\cos(k_{n}))\left(a^{\dagger}_{k_{n}}a_{k_{n}}-a_{-k_{n}}a^{\dagger}_{-k_{n}}\right)+\Delta\left(e^{-ik_{n}}a^{\dagger}_{-k_{n}}a^{\dagger}_{k_{n}}+e^{ik_{n}}a_{k_{n}}a_{-k_{n}}\right)\right].\end{split}$
(13)

The values of $k_{n}$ are fixed by requiring that the momentum-space creation
and annihilation operators satisfy the correct anti-commutation relations as
well as imposing the desired boundary condition on $c_{j}$. The generalized
boundary condition,

$c_{L+1}=e^{i\phi}c_{1},$ (14)

with $e^{i\phi}=1\textnormal{ and }e^{i\phi}=-1$ denoting periodic and anti-
periodic boundary conditions respectively, corresponds to the set of momenta

$k_{n}=\frac{2\pi}{L}\left(-\left\lceil\frac{L-1}{2}\right\rceil+n+\frac{\phi}{2\pi}\right)\
\ \ ;\ \ \ n=0,1,2,\cdots L-1\,.$ (15)

For even $L$, the Kitaev chain with the anti-periodic boundary condition
(APBC) may be written as

$H=-\sum_{k_{n}>0}\left[2(\mu+J\cos(k_{n}))J_{0}^{(k_{n})}-i\Delta\sin(k_{n})\left(J_{+}^{(k_{n})}-J_{-}^{(k_{n})}\right)\right].$
(16)

Closely related expressions for periodic boundary conditions and for odd $L$
can be found in Appendix A of Kawabata:2017zsb . In this form of the
Hamiltonian we defined the ladder operators

$\displaystyle J_{0}^{(k_{n})}$ $\displaystyle\equiv$
$\displaystyle\frac{1}{2}\left(a^{\dagger}_{k_{n}}a_{k_{n}}-a_{-k_{n}}a^{\dagger}_{-k_{n}}\right)\,,$
$\displaystyle J_{+}^{(k_{n})}$ $\displaystyle\equiv$ $\displaystyle
a^{\dagger}_{k_{n}}a^{\dagger}_{-k_{n}}\,,$ (17) $\displaystyle
J_{-}^{(k_{n})}$ $\displaystyle\equiv$ $\displaystyle a_{-k_{n}}a_{k_{n}}\,,$

that satisfy the $su(2)$ algebra

$\displaystyle\left[J_{0}^{(k_{n})},J_{\pm}^{(k_{n^{\prime}})}\right]$
$\displaystyle=$ $\displaystyle\pm\delta_{n,n^{\prime}}J_{\pm}^{(k_{n})}\,,$
$\displaystyle\left[J_{+}^{(k_{n})},J_{-}^{(k_{n^{\prime}})}\right]$
$\displaystyle=$ $\displaystyle 2\delta_{n,n^{\prime}}J_{0}^{(k_{n})}\,.$

A key point to note is that, after defining,

$\vec{R}(k_{n})\equiv\begin{pmatrix}R_{1}\\\ R_{2}\\\ R_{3}\\\
\end{pmatrix}=\begin{pmatrix}0\\\ \Delta\sin k_{n}\\\ \mu+J\cos k_{n}\\\
\end{pmatrix},$ (18)

the Hamiltonian (16) is essentially of the same form

$H=-\sum_{k_{n}>0}\left[2R_{3}J_{0}^{(k_{n})}-iR_{2}(J_{+}^{(k_{n})}-J_{-}^{(k_{n})})\right],$
(19)

as that studied in Caputa:2022eye , with the only differences being the values
of the coefficients $R_{2},R_{3}$ as well as only summing over positive values
of $k$. The latter feature may be anticipated from the map between the SSH
Hamiltonian and the periodic Kitaev chain Bandyopadhyay:2021dei .

The full excitation spectrum of the Hamiltonian may be represented in terms of
$SU(2)$ coherent states444For each $k$ the states
$|\psi_{+}\rangle=\cos(\frac{\varphi_{k}}{2})e^{-i\tan(\frac{\varphi_{k}}{2})J_{+}^{(k)}}|\frac{1}{2},-\frac{1}{2}\rangle$
and
$|\psi_{-}\rangle=\sin(\frac{\varphi_{k}}{2})e^{i\cot(\frac{\varphi_{k}}{2})J_{+}^{(k)}}|\frac{1}{2},-\frac{1}{2}\rangle$
are eigenstates with eigenvalues
$\epsilon_{\pm}=\pm\sqrt{(\mu+J\cos\frac{\varphi_{k}}{2})^{2}+\left(\Delta\sin\frac{\varphi_{k}}{2}\right)^{2}}$.
and, in particular, the ground state is given by

$|\Omega\rangle=\prod_{k_{n}>0}\sin\left|\frac{\varphi_{k_{n}}}{2}\right|\exp\left[e^{-i\psi_{k_{n}}}\cot\frac{\varphi_{k_{n}}}{2}J_{+}^{(k_{n})}\right]\left|\frac{1}{2},-\frac{1}{2}\right\rangle_{k_{n}}\,,$
(20)

where we have introduced the spherical coordinates
$(R(k_{n}),\varphi_{k_{n}},\psi_{k_{n}})$ to rewrite the parameter vector
$\vec{R}(k_{n})$. In particular,

$\left\\{\begin{matrix}R(k_{n})=\left|\vec{R}(k_{n})\right|\\\
\varphi_{k_{n}}=\arctan\frac{|\Delta|\sin k_{n}}{\mu+J\cos k_{n}}\\\
\psi_{k_{n}}=\frac{\pi}{2}-\operatorname{arg}\Delta\\\
\end{matrix}\right.,\qquad k_{n}>0\,.$ (21)

Note that these relations remain true for complex $\Delta$ as well. To compare
our conventions with Liu:2019aji , set

$\theta_{k_{n}}=\operatorname{sgn}\Delta\cdot\frac{\pi-\varphi_{k_{n}}}{2}\,.$

The reference state above is a tensor product of the lowest weight $su(2)$
states in the spin-$\frac{1}{2}$ representation, for each positive value of
$k_{n}$.

Before proceeding, a few points are in order.

  * •

The first is a comment on some of the existing literature on this topic. The
periodic Kitaev chain with even $L$ and both short- and long-range pairing has
been studied in Liu:2019aji  where it was determined that Nielsen’s approach
to circuit complexity may be used to detect the topological phase transition.
In this paper we will instead study the spread complexity to determine whether
it is also sensitive to this transition. Another model exhibiting a
topological phase transition, the SSH model Su:1979ua , has been studied using
both circuit complexity Ali:2018aon  and spread complexity Caputa:2022eye .
Taken together, these results seem to point to spread complexity being
sensitive to the topological phase transition, while Nielsen complexity is
not. This would appear to be in contrast to our observations of the
topological transition in the Kitaev chain. There are a few possible
explanations of this, the simplest of which is that in Ali:2018aon  the
amplitude of the fluctuations in the early time regime, which is essential to
detect the transition as argued in Caputa:2022eye , was not probed carefully
enough.

  * •

The second concerns a generalization of the Kitaev chain to account for
twisted boundary conditions. The associated Hamiltonian and site-dependent
chemical potential are Kawabata:2017zsb

$\displaystyle H$ $\displaystyle=$
$\displaystyle\sum_{j=1}^{L-1}\left[-\frac{J}{2}\left(c_{j}^{\dagger}c_{j+1}+c_{j+1}^{\dagger}c_{j}\right)+\frac{\Delta}{2}\left(c_{j}c_{j+1}+c_{j+1}^{\dagger}c_{j}^{\dagger}\right)\right]-\sum_{j=1}^{L}\mu_{j}\left(c_{j}^{\dagger}c_{j}-\frac{1}{2}\right)$
$\displaystyle+b\left[-\frac{J}{2}(e^{i\phi_{1}}c_{L}^{\dagger}c_{1}+e^{-i\phi_{1}}c_{1}^{\dagger}c_{L})+\frac{\Delta}{2}(e^{i\phi_{2}}c_{L}c_{1}+e^{-i\phi_{2}}c_{1}^{\dagger}c_{L}^{\dagger})\right],$
$\displaystyle\mu_{j}$ $\displaystyle=$
$\displaystyle\mu+\mu(a-1)(\delta_{j,1}+\delta_{j,L}).$ (22)

The bulk chemical potential, hopping amplitude and superconducting pairing
strength are as before, while the chemical potential at the boundary sites
$i=1,L$ has a modified form $\mu_{1}=\mu_{L}=a\mu$. In addition, the terms
connecting sites $1$ and $L$ pick up additional phases $\phi_{1},\phi_{2}$.
These may run over an interval of $2\pi$ and we will find it convenient to
parameterise the interval as $\phi_{i}\in[-\pi,\pi)$.

As usual, the model exhibits a topological phase transition when the ground
state energy becomes degenerate. As shown in Kawabata:2017zsb , this coincides
with the existence of Majorana zero modes. Such zero modes are not present for
arbitrary choices of the system parameters. This model should thus allow one
to assess spread complexity as a probe of the topological phase transition.
Essentially, one would want a probe to diagnose the phase transition if and
only if it is present. To perform the computation one would have to obtain,
numerically or otherwise, the Krylov basis. We leave this detailed study to
future work but provide some promising approximate results in this direction
in Appendix B.

##  4 Spread complexity as a probe of topological phase transitions

As previously discussed, eigenstates of the Kitaev Hamiltonian (19) may be
fully described in terms of $SU(2)$ coherent states. In addition, spread
complexity is a function of the Hamiltonian parameters and the choice of
reference and target state which determine the circuit. In what follows, we
consider circuits connecting the ground state of the Kitaev chain (taken as
the target state) to a number of natural and universally definable choices of
reference states. Following Caputa:2022eye  we compute the spread complexity
for each of these circuits and show its sensitivity to the topological phase
transition in all cases. In this sense, this demonstrates the robustness of
spread complexity as a probe of topological phase transitions, at least for
the Kitaev chain and the models related to it. Unless otherwise stated we work
with the Hamiltonian (16) valid for even $L$ and anti-periodic boundary
conditions.

###  4.1 Circuit 1

All the circuits considered here and below will take the Kitaev chain ground
state as the target state and we distinguish between the different circuits by
their choice of reference state. For convenience, let us rewrite the ground
state (20) utilising the Baker-Campbell-Hausdorff identity, as

$\displaystyle|\Omega(s=1)\rangle$ $\displaystyle=$
$\displaystyle\bigotimes_{k>0}|\Omega_{k}(s=1)\rangle,$ (23)
$\displaystyle|\Omega_{k}(s=1)\rangle$ $\displaystyle=$
$\displaystyle\exp\left[\frac{\pi-\varphi_{k}}{2}e^{-i\left(\frac{\pi}{2}-\operatorname{arg}\Delta\right)}J_{+}^{(k)}-\frac{\pi-\varphi_{k}}{2}e^{i\left(\frac{\pi}{2}-\operatorname{arg}\Delta\right)}J_{-}^{(k)}\right]\left|\frac{1}{2},-\frac{1}{2}\right\rangle_{k},$
(24) $\displaystyle=$
$\displaystyle\exp\left[-i\frac{\pi-\varphi_{k}}{2}\left(e^{i\operatorname{arg}\Delta}J_{+}^{(k)}+e^{-i\operatorname{arg}\Delta}J_{-}^{(k)}\right)\right]\left|\frac{1}{2},-\frac{1}{2}\right\rangle_{k}.$

As a first example of reference state, we will take

$|\Omega_{k}(s=0)\rangle=\left|\frac{1}{2},\ \operatorname{sgn}\left(\mu+J\cos
k\right)\cdot\frac{1}{2}\right\rangle_{k}.$ (25)

In a sense, this reference state is a “natural” one to choose, since it is
nothing but the ground state of the free fermionic Hamiltonian (or the mass
part of the Hamiltonian) in (16), i.e., it is the ground state of

$H_{M}=-\sum_{k>0}2R_{3}J_{0}^{(k)}=-\sum_{k>0}2|R_{3}|\left(\operatorname{sgn}R_{3}\cdot
J_{0}^{(k)}\right),$ (26)

with $R_{3}=\mu+J\cos k$ and where in the last expression above, we have
expressed the Hamiltonian in a form that makes the ground state more obvious.
Note that (25) is equivalent to

$|\Omega_{k}(s=0)\rangle=\exp\left[-i\frac{\pi}{2}\Theta\left(\mu+J\cos
k\right)\cdot\left(J_{+}^{(k)}+J_{-}^{(k)}\right)\right]\left|\frac{1}{2},-\frac{1}{2}\right\rangle_{k},$
(27)

where we have introduced the step function $\Theta$.

The first circuit we consider, connecting the free fermion ground state to the
Kitaev chain ground state, thus takes the form

$|\Omega_{k}(s)\rangle=\exp\left[-is\left(\frac{\pi-\varphi_{k}}{2}\operatorname{sgn}\Delta-\frac{\pi}{2}\Theta\left(\mu+J\cos
k\right)\right)\left(J_{+}^{(k)}+J_{-}^{(k)}\right)\right]|\Omega_{k}(s=0)\rangle.$
(28)

The spread complexity is then given by

$\displaystyle\mathcal{C}_{k}(s=1)$ $\displaystyle=$
$\displaystyle\sin^{2}\left(\theta_{k}-\frac{\pi}{2}\Theta\left(\mu+J\cos
k\right)\right)=\frac{1}{2}\left(1+\operatorname{sgn}\left(\mu+J\cos
k\right)\cdot\cos 2\theta_{k}\right)$ (29) $\displaystyle=$
$\displaystyle\frac{1}{2}\left(1-\operatorname{sgn}\left(\mu+J\cos
k\right)\cdot\cos\varphi_{k}\right)$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left(1-\frac{|\mu+J\cos(k)|}{\sqrt{(\mu+J\cos(k))^{2}+\Delta^{2}\sin^{2}(k)}}\right).$

To find the total spread complexity we have to sum over the discrete values of
$k_{n}$ for finite $L$, which becomes an integral in the continuum limit

$C(J,\mu,\Delta)=\frac{1}{L}\sum_{n>0}C_{k_{n}}(s=1)\rightarrow\frac{1}{\pi}\int_{0}^{\pi}dk\,C_{k}(s=1)\,.$
(30)

Figure 1: The spread complexity in the continuum limit for the circuit
connecting the free fermion ground state to the Kitaev chain ground state. We
have chosen $J=1$. When $|\mu|<1$ the system is in the topological phase and
the spread complexity is a $\Delta$-dependent constant.

Note that the normalization is with respect to the total number of sites in
the system. Although we have not found an analytic expression for the
continuum integral, it is not difficult to see that

$C(J,|\mu|<|J|,\Delta)=\frac{1}{2}-\frac{1}{\pi}\int_{0}^{\frac{\pi}{2}}\sin(\varphi_{k})\tan^{-1}\left(\frac{\tan(\varphi_{k})}{|\Delta|}\right)d\varphi_{k}\,,$
(31)

which is $\mu$-independent. In the continuum limit we thus find,
qualitatively, the same result as Caputa:2022eye  where the spread complexity
of the circuit connecting the free fermion ground state and Kitaev chain
ground state is constant in the topological phase. Consequently, we conclude
that spread complexity is indeed sensitive to the phase transition, as can be
seen in figure (1).

To compare with Caputa:2022eye , we also plot the spread complexity for a few
values of $\Delta$ in Fig. (2(a)). The spread complexity takes a constant
$\Delta$-dependent value in the topological phase. At finite $L$ the spread
complexity is not constant in the topological phase, but rather resembles a
gluing of a discrete number of arc segments as $\mu$ is changed within this
phase, as in figure (2(b)). As $L$ is increased, the number of these segments
increase and the profile approximates the continuum limit constant value that
characterizes the topological phase.

(a)

(b)

Figure 2: a) The spread complexity for the circuit connecting the free fermion
ground state and Kitaev chain ground state as a function of the chemical
potential, with $J=1$, and in the continuum limit. The complexity is constant
in the region $|\mu|<1$. b) The spread complexity for the circuit connecting
the free fermion ground state and Kitaev chain ground state as a function of
the chemical potential, with $J=1$ and $\Delta=2$. The spread complexity has
been divided by $\sqrt{L}$ for display purposes. The top curve corresponds to
$L=40$, the middle $L=20$ and the bottom $L=10$.

Another feature we report on is the behavior of the spread complexity as
$\Delta$ approaches zero. From (29), it is clear that the spread complexity is
both continuous and vanishes at $\Delta=0$. The behavior of the first
derivative $\frac{\partial C(J,\mu,\Delta)}{\partial\Delta}$, however, is
different in the topological and trivial phases as $\Delta$ approaches zero.
Specifically, when $|\mu|>|J|$ the first derivative approaches zero as
$\Delta$ approaches zero while for $|\mu|<|J|$ the first derivative approaches
a constant value. Viewed as a function of $\Delta$, running over positive and
negative values, the spread complexity develops a cusp in the topological
phase at $\Delta=0$. While we have no physical explanation for this behaviour
at present, it is reminiscent of the behaviour of the order parameter of a
conventional second order phase transition. In figure (3) we plot the behavior
of the derivative for two values of $\mu$ near the phase transition to
illustrate this observation.

Figure 3: The derivative of spread complexity with respect to $\Delta$ in the
continuum limit for the circuit connecting the free fermion ground state and
Kitaev chain ground state. We have set $J=1$. The blue line has $\mu=0.98$,
the orange $\mu=1.02$ and the green $\mu=1.1$. When crossing the topological
phase transition points at $|\mu|=1$ the derivative develops a discontinuity.

###  4.2 Circuit 2

The spread complexity of the circuit in section 4.1 exhibits a plateau-like
feature when the system is in the topological phase. This is reminiscent of
results obtained for the SSH model Caputa:2022eye , where the spread
complexity also plateaued in the topological phase. An undesirable feature of
the first circuit, however, is that the reference state (25) is
$\mu$-dependent. A natural question to ask is whether we can make a
$(\mu,\Delta)$-independent choice for the reference state such that the spread
complexity exhibits a plateau in topological phase? The answer it turns out is
yes, if we select reference state to be

$|\Omega_{k}(s=0)\rangle=\exp\left[\frac{\pi}{4}e^{-i\frac{\pi}{2}}J_{+}^{(k)}-\frac{\pi}{4}e^{+i\frac{\pi}{2}}J_{-}^{(k)}\right]\left|\frac{1}{2},-\frac{1}{2}\right\rangle_{k}.$
(32)

This state can be interpreted physically as the ground state of the Cooper-
pair (fermion interaction) part of the Kitaev Hamiltonian. With this choice of
reference state the circuit takes the form

$|\Omega_{k}(s)\rangle=\exp\left[-i\cdot
s\cdot\left(\operatorname{sgn}\Delta\cdot\frac{\pi-\varphi_{k}}{2}-\frac{\pi}{4}\right)\left(J_{+}^{(k)}+J_{-}^{(k)}\right)\right]|\Omega_{k}(s=0)\rangle\,,$
(33)

from which we compute the spread complexity to be

$\mathcal{C}_{k}(s=1)=\sin^{2}\left(\operatorname{sgn}\Delta\cdot\frac{\pi-\varphi_{k}}{2}-\frac{\pi}{4}\right)=\frac{1}{2}\left(1-\frac{\Delta\sin
k}{\sqrt{(\mu+J\cos(k))^{2}+\Delta^{2}\sin^{2}(k)}}\right)\,.$ (34)

Integrating by parts, we find that in the regime $|\mu|<|J|$ the total spread
complexity is given by

$\mathcal{C}(J,|\mu|<|J|,\Delta)=\frac{1-\operatorname{sgn}\Delta}{2}+\frac{1}{\pi}\int_{0}^{\pi/2}\,d\varphi_{k}\,\cos\varphi_{k}\,\tan^{-1}\frac{\tan\varphi_{k}}{\Delta}\,.$
(35)

In the topological phase the spread complexity is therefore independent of the
chemical potential, $\mu$. The integral can be performed analytically to give

$\mathcal{C}(J,|\mu|<|J|,\Delta)=\frac{1}{2}-\frac{1}{\pi}\frac{\Delta\tan^{-1}\sqrt{\Delta^{2}-1}}{\sqrt{\Delta^{2}-1}}\,.$
(36)

We plot the spread complexity as a function of $\mu$ for various choices of
$\Delta$ in Fig. (4). The plateau feature of spread complexity in the
topological phase can be clearly observed. In fig. (5) we plot the dependence
on $\Delta$ of spread complexity inside the topological phase i.e. the height
of the plateau. The red dots on the figure coincide with the plateau value for
the SSH model Caputa:2022eye , at $\Delta=-1$. At these points the numerator
and denominator of (36) tend to zero, so that the value is obtained by taking
the limit $\Delta^{2}\rightarrow 1$.

As in the previous case, the derivatives of spread complexity, in the
topological phase, also exhibit noteworthy behaviour. Specifically, the
derivative of complexity with respect to $\Delta$,

$\pi\frac{d}{d\Delta}\mathcal{C}(J,|\mu|<|J|,\Delta)=\frac{1}{1-\Delta^{2}}+\frac{\tan^{-1}\left(\sqrt{\Delta^{2}-1}\right)}{\left(\Delta^{2}-1\right)^{3/2}},$
(37)

is continuous everywhere, including $|\Delta|=1$. The second derivative is
also continuous at $|\Delta|=1$ but diverges at $\Delta=0$, since

$\pi\frac{d^{2}}{d\Delta^{2}}\mathcal{C}(J,|\mu|<|J|,\Delta)=\frac{2\Delta^{2}+1}{\Delta\left(\Delta^{2}-1\right)^{2}}-\frac{3\Delta\tan^{-1}\left(\sqrt{\Delta^{2}-1}\right)}{\left(\Delta^{2}-1\right)^{5/2}}\sim\Delta^{-1}\,,$
(38)

as $\Delta\to 0$. Again, note that this curious divergence of the second
derivative is reminiscent of a third order (conventional) phase transition and
deserves further investigation that we leave for future explorations.

Figure 4: Complexity $\mathcal{C}(s=1;\mu,\Delta)$ as a function of $\mu$ for
various choices of $\Delta$ ($\Delta=-2,-1,-1/2,1/2,1,2$) for the circuit
connecting the Cooper-pair ground state and Kitaev chain ground state. Between
the two vertical gridlines $\mu=\pm 1$, spread complexity is
$\mu$-independent. The two dashed horizontal gridlines are the analytical
results of $\mathcal{C}(|\mu|<1,\Delta=\pm 1)$, respectively. For
$|\mu|\to\infty$, the spread complexities of various $\Delta$s approach the
dotted horizontal gridline $\mathcal{C}=1/2$. Figure 5: Complexity
$\mathcal{C}(s=1;|\mu|<1,\Delta)$ as a function of $\Delta$ for the circuit
connecting the Cooper-pair ground state and Kitaev chain ground state. The two
red dots are the analytical results of $\mathcal{C}(|\mu|<1,\Delta=\pm 1)$,
respectively.

###  4.3 Circuit 3

Finally, let’s consider the circuit connecting the fermion vacuum to the
Kitaev chain ground state. The fermion vacuum, with its complete absence of
particle states, is in some sense the most straightforward choice of reference
state. In the $SU(2)$ representation it is the tensor product (over momenta)
of lowest weight spin-$\frac{1}{2}$ states. For each value of $k$ the circuit
takes the form

$|\Omega_{k}(s)\rangle=\exp\left[s\cdot\left(\frac{\pi-\varphi_{k}}{2}e^{-i\left(\frac{\pi}{2}-\operatorname{arg}\Delta\right)}J_{+}^{(k)}-\frac{\pi-\varphi_{k}}{2}e^{i\left(\frac{\pi}{2}-\operatorname{arg}\Delta\right)}J_{-}^{(k)}\right)\right]\left|\frac{1}{2},-\frac{1}{2}\right\rangle_{k},$
(39)

where we have again made use of Baker-Campbell-Hausdorff identities to rewrite
the ground state (20). We obtain the spread complexity as

$\mathcal{C}_{k}(s=1)=\sin^{2}\frac{\pi-\varphi_{k}}{2}=\frac{1}{2}\left(1+\frac{\mu+J\cos(k)}{\sqrt{(\mu+J\cos(k))^{2}+\Delta^{2}\sin^{2}(k)}}\right).$
(40)

The derivative of the spread complexity with respect to $\mu$,

$\frac{\partial}{\partial\mu}\mathcal{C}_{k}(s=1)=\frac{\Delta^{2}\sin^{2}(k)}{2((\mu+J\cos(k))^{2}+\Delta^{2}\sin^{2}(k))^{\frac{3}{2}}}\,,$
(41)

will also prove important. As before, we need to sum over the values of $k$ in
order to obtain the spread complexity for this circuit. In figure (6(a)), we
plot the spread complexity as a function of $\mu$ in the continuum limit.
Surprisingly, unlike for the previous circuits, the complexity monotonically
increases with $\mu$. In particular we do not observe any plateau in the
spread complexity in the topological phase. The spread complexity still has a
noteworthy feature at the phase transition point. Specifically, the first
derivative of complexity diverges at the phase transition points, as in figure
(6(b)). The divergence at $\mu=-1$ is due to low momentum contributions while
the divergence at $\mu=1$ is due to the high momentum contributions,

$\displaystyle\left.\frac{\partial}{\partial\mu}\mathcal{C}_{k}(s=1)\right|_{\mu=1}$
$\displaystyle=$ $\displaystyle\frac{1}{2|\Delta|k}+\mathcal{O}(k)\,,$
$\displaystyle\left.\frac{\partial}{\partial\mu}\mathcal{C}_{k}(s=1)\right|_{\mu=-1}$
$\displaystyle=$
$\displaystyle\frac{1}{2|\Delta|(k-\pi)}+\mathcal{O}(k-\pi)\,.$ (42)

For the intermediate regime $|\mu|<1$ there always exists a pole at some value
of $k\in\left[0,\pi\right]$ at which the integral diverges. Evidently, even
though we do not observe a plataeu in the topological phase, the spread
complexity for this circuit is able to distinguish the topological phase
transition in the chain.

(a)

(b)

Figure 6: a) The spread complexity for the circuit connecting the vacuum and
the Kitaev chain ground state in the continuum limit. We have set $J=1$. The
spread complexity is monotonically increasing as a function of $\mu$. b) The
inverse of the derivative of the spread complexity as a function of $\mu$ for
$\Delta=\frac{1}{2}$ and $J=1$. At the phase transition the derivative
diverges so that its inverse approaches zero.

At finite $L$ we observe that the spread complexity has the same schematic
form depicted in (6(a)). The derivative, however, remains finite, with a sharp
change in the sign of the second derivative. As $L$ increases the derivative
increases, reaching the diverging value in the $L\rightarrow\infty$ limit. In
this circuit the behavior around $\Delta=0$ is also an alternate diagnostic of
the topological phase transition. Using (40) we may compute

$\displaystyle\left.\frac{\partial^{2}}{\partial\Delta^{2}}C_{k}(s=1)\right|_{\mu=-1,\Delta=0}$
$\displaystyle=$ $\displaystyle\frac{1}{2}\cot^{2}\left(\frac{k}{2}\right),$
$\displaystyle\left.\frac{\partial^{2}}{\partial\Delta^{2}}C_{k}(s=1)\right|_{\mu=-1,\Delta=0}$
$\displaystyle=$ $\displaystyle-\frac{1}{2}\tan^{2}\left(\frac{k}{2}\right),$
(43)

which diverges when integrated over the interval $k\in\left[0,\pi\right]$. In
the region $|\mu|>1$ we find that the second derivative remains finite while
in the region $|\mu|\leq 1$ the integral diverges. As in the first circuit, a
pole develops in the second derivative at some value of $k$ inside the
integration range, leading to a diverging integral.

###  4.4 Odd $L$ and Periodic boundary conditions

The results in the previous three subsections were obtained by studying the
APBC Kitaev chain with an even number of sites and afterwards taking the limit
to an infinite number of sites. One may wonder how this choice of boundary
condition and the parity of sites affect the results. In this subsection we
briefly explain why the large $L$ limit results are insensitive to these
choices.

In the case of odd $L$ and/or a periodic boundary condition on the Kitaev
chain, the momenta (15) can include $k=0$ and/or $k=\pi$. Since, for these
values, $a_{k}=a_{-k}$, this situation requires special attention. This adds
one of or both of the terms

$\displaystyle H_{0}$ $\displaystyle=$
$\displaystyle-(\mu+J)(a_{0}^{\dagger}a_{0}-a_{0}a_{0}^{\dagger}),$
$\displaystyle H_{\pi}$ $\displaystyle=$
$\displaystyle-(\mu-J)(a_{\pi}^{\dagger}a_{\pi}-a_{\pi}a_{\pi}^{\dagger}),$
(44)

to the Hamiltonian (16), see e.g. Appendix A of Kawabata:2017zsb . The
eigenstates of $H_{0},H_{\pi}$ is simply the number operator basis i.e. a
single fermion or no fermion. The ground state for $H_{0}$ and $H_{\pi}$ is
given by

$\displaystyle|\Omega\rangle_{0}$ $\displaystyle=$
$\displaystyle|0\rangle_{0}+\Theta(\mu+J)\left(a_{0}^{\dagger}|0\rangle_{0}-|0\rangle_{0}\right),$
$\displaystyle|\Omega\rangle_{\pi}$ $\displaystyle=$
$\displaystyle|0\rangle_{\pi}+\Theta(\mu-J)\left(a_{\pi}^{\dagger}|0\rangle_{\pi}-|0\rangle_{\pi}\right),$
(45)

respectively. The spread complexity can be computed using the Heisenberg
algebra formulae of Caputa:2021sib  which adds a factor of $1$ to the total
complexity if the ground state differs from the reference state by a single
fermion.

In the large $L$ limit this contribution is subleading and the results of the
previous subsections hold without modification. At finite $L$ we need to add
$1$ to the complexity at $\mu=-J$ ($\mu=J$) if $H_{0}$ ($H_{\pi}$) is included
for the boundary condition under consideration. This shift has the effect of
transforming the finite-$L$ spread complexity from a continuous function to a
discontinuous one at these points.

##  5 Discussion

Computational complexity, a well-known entity in the context of computer
science and associated mathematical subfields like graph theory, is still very
much a shiny new concept in high energy physics. While the primary motivation
for its introduction in the latter was the understanding of the black hole
information-loss paradox, it has proven far more versatile with recent
applications to the emergence of spacetime in quantum gravity, early universe
cosmology Bhattacharyya:2020rpy ; Bhattacharyya:2020kgu ; Haque:2021kdm ;
Haque:2021hyw  and even condensed matter physics. It is this last aspect that
is the focus of this article.

Specifically, this paper details our study of a particular variant of quantum
complexity called spread complexity, introduced in Balasubramanian:2022tpr ,
applied to Kitaev’s 1-dimensional p-wave superconductor. The Kitaev chain is a
prototypical example of a system exhibiting topological phase transitions, and
an important laboratory for tests of novel features of quantum systems such as
topological order and quantum phases of matter. Building on the methodology of
Caputa:2022eye , we show here that spread complexity, based on the associated
Krylov basis of the Kitaev chain, is sensitive to the topological phase
transition exhibited by the model. Our results are in concordance with, and
provides further support for, those of Caputa:2022eye , where results for the
SSH model were also in the affirmative for spread complexity as a novel probe
of the topological phase transition.

Furthermore we have demonstrated a level of robustness of these results by
computing the spread complexity for circuits of differing reference states
connected to the Kitaev ground state taken as the fixed target state. When
using the $\Delta=0$ or $J=0$ ground state as reference state we have shown
that, in the topological phase, the spread complexity is a constant as a
function of $\mu$. Additionally, its second derivative with respect to
$\Delta$ exhibits divergences when the system is in the topological phase.
When using the vacuum as reference state the spread complexity does not
exhibit the plateau feature, but its first derivative with respect to $\mu$
diverges at the transition point. In all these cases the spread complexity
thus does, in some way, capture the transition from the trivial to topological
phase. In essence, the statement that we are making is that two states
realizing distinct topological phases, certainly in the Kitaev chain and the
SSH models cannot be mapped to one another via low-depth circuits.

The circuits that we have considered connect simple choices of reference state
to the ground state of the Kitaev chain. Even though the spread complexity
does indeed capture the transition from trivial phase to topological phase in
some way for each of our circuit choices, it is clear that not all choices
exhibit a plateau in the topological phase. Based on results for the SSH model
Caputa:2022eye  this may seem surprising, and it leads to a natural question:
which circuit choices give rise to a spread complexity plateau? Presumably, an
understanding of this feature should be related to the symmetries of states in
different phases and it may be possible to classify circuits according to
whether they preserve or break these symmetries. Hopefully, elucidating this
symmetry classification may provide some hint as to the general choices of
circuits that would display the plateau feature in systems exhibiting
topological phase transitions.

The utility of spread complexity as an efficient probe of quantum phase
transitions can be extended and explored in several different directions:

  * •

Given the divergence of fluctuations at a quantum critical point, it might
intuitively be expected that the complexity of say the ground state of a
system undergoing a quantum phase transition diverges at the transition. This
article reflects our attempts to sharpen this intuition and, like its
predecessors in the literature Ali:2018aon ; Liu:2019aji ; Caputa:2021sib
focused on a particular type of quantum phase transition in which the
transition is between states realizing distinct topological phases. Here we
can indeed make a crisp statement about the complexity (and its derivatives)
at the transition. However, topological phase transitions are just one in an
ever-expanding list of quantum phase transitions that include conventional
(Landau) order, Berezinskii-Kosterlitz-Thouless-transitions Kosterlitz_1973
and the more exotic deconfined critical points of Senthil_2004 . It would be
extremely exciting if spread complexity, or the closely related K-complexity
offered a window into these phase transitions as well. Indeed one tantalising
possibility, which we leave for future work, might lie in the well-known
relation between the Kitaev chain and the transverse field Ising model that
map the topological order of the former to conventional order in the latter.

  * •

The computations in this article rely strongly on us being able to represent
the Kitaev Hamiltonian in a basis of $su(2)$ generators. But this is not
always possible; a fact that is demonstrable even in the Kitaev chain itself.
For example, when the (anti)periodic boundary conditions studied above are
extended to more general twisted boundary conditions, the Hamiltonian can no
longer be written in terms of $su(2)$ generators. Understanding how to compute
the spread complexity in such situations would be an important development of
this computational framework. As a first step in this direction, in Appendix
B, as a preview of forthcoming work, we give some preliminary results for the
Kitaev chain with twisted boundary conditions.

###### Acknowledgements.

JM would like to acknowledge support from the ICTP through the Associates
Programme and from the Simons Foundation through grant number 284558FY19 as
well as funding from the South African Research Chairs Initiative funded by
the National Research Foundation and Department of Science and Technology. NG
is supported by a Faculty of Science PhD Fellowship. H.J.R.vZ is supported by
the “Quantum Technologies for Sustainable Devlopment" grant from the National
Institute for Theoretical and Computational Sciences (NITHECS). Work of PC and
SL is supported by “Polish Returns 2019” grant of the National Agency for
Academic Exchange (NAWA) PPN/PPO/2019/1/00010/U/0001 and Sonata Bis 9
2019/34/E/ST2/00123 grant from the National Science Centre, NCN.

## References

  * (1) A. Y. Kitaev, Unpaired majorana fermions in quantum wires, Physics-Uspekhi 44 (oct, 2001) 131–136.
  * (2) M. Leijnse and K. Flensberg, Introduction to topological superconductivity and Majorana fermions, Semicond. Sci. Tech. 27 (2012) 124003, [arXiv:1206.1736].
  * (3) M. Greiter, V. Schnells, and R. Thomale, The 1d ising model and topological order in the kitaev chain, Annals of Physics 351 (02, 2014) [arXiv:1402.5262].
  * (4) W. P. Su, J. R. Schrieffer, and A. J. Heeger, Solitons in polyacetylene, Phys. Rev. Lett. 42 (1979) 1698–1701.
  * (5) F. Franchini, An introduction to integrable techniques for one-dimensional quantum systems, vol. 940.  Springer, 2017.
  * (6) M. A. Nielsen and I. L. Chuang, Quantum Computation and Quantum Information.  Cambridge University Press, 2000.
  * (7) J. Watrous, Quantum computational complexity, Computational Complexity: Theory, Techniques, and Applications (05, 2008) [arXiv:0804.3401].
  * (8) J. M. Maldacena, The Large N limit of superconformal field theories and supergravity, Adv. Theor. Math. Phys. 2 (1998) 231–252, [hep-th/9711200].
  * (9) A. R. Brown, D. A. Roberts, L. Susskind, B. Swingle, and Y. Zhao, Holographic Complexity Equals Bulk Action?, Phys. Rev. Lett. 116 (2016), no. 19 191301, [arXiv:1509.07876].
  * (10) N. Y. Halpern, N. B. T. Kothakonda, J. Haferkamp, A. Munson, J. Eisert, and P. Faist, Resource theory of quantum uncomplexity, arXiv:2110.11371.
  * (11) D. Harlow and P. Hayden, Quantum Computation vs. Firewalls, JHEP 06 (2013) 085, [arXiv:1301.4504].
  * (12) S. Aaronson, The Complexity of Quantum States and Transformations: From Quantum Money to Black Holes, 7, 2016.  arXiv:1607.05256.
  * (13) L. Susskind, Computational Complexity and Black Hole Horizons, Fortsch. Phys. 64 (2016) 24–43, [arXiv:1403.5695]. [Addendum: Fortsch.Phys. 64, 44–48 (2016)].
  * (14) D. Stanford and L. Susskind, Complexity and Shock Wave Geometries, Phys. Rev. D 90 (2014), no. 12 126007, [arXiv:1406.2678].
  * (15) J. Maldacena, S. H. Shenker, and D. Stanford, A bound on chaos, JHEP 08 (2016) 106, [arXiv:1503.01409].
  * (16) A. Bouland, B. Fefferman, and U. Vazirani, Computational pseudorandomness, the wormhole growth paradox, and constraints on the AdS/CFT duality, arXiv:1910.14646.
  * (17) S. Chapman and G. Policastro, Quantum computational complexity from quantum information to black holes and back, Eur. Phys. J. C 82 (2022), no. 2 128, [arXiv:2110.14672].
  * (18) B. Chen, B. Czech, and Z.-z. Wang, Quantum information in holographic duality, Rept. Prog. Phys. 85 (2022), no. 4 046001, [arXiv:2108.09188].
  * (19) M. A. Nielsen, A geometric approach to quantum circuit lower bounds, 2005\.
  * (20) M. A. Nielsen, M. R. Dowling, M. Gu, and A. C. Doherty, Quantum computation as geometry, Science 311 (feb, 2006) 1133–1135, [quant-ph/0603161].
  * (21) M. R. Dowling and M. A. Nielsen, The geometry of quantum computation, 2007\.
  * (22) T. Ali, A. Bhattacharyya, S. Shajidul Haque, E. H. Kim, and N. Moynihan, Post-Quench Evolution of Complexity and Entanglement in a Topological System, Phys. Lett. B 811 (2020) 135919, [arXiv:1811.05985].
  * (23) F. Liu, S. Whitsitt, J. B. Curtis, R. Lundgren, P. Titum, Z.-C. Yang, J. R. Garrison, and A. V. Gorshkov, Circuit complexity across a topological phase transition, Phys. Rev. Res. 2 (2020), no. 1 013323, [arXiv:1902.10720].
  * (24) Z. Xiong, D.-X. Yao, and Z. Yan, Nonanalyticity of circuit complexity across topological phase transitions, Phys. Rev. B 101 (2020), no. 17 174305, [arXiv:1906.11279].
  * (25) P. Caputa and S. Liu, Quantum complexity and topological phases of matter, arXiv:2205.05688.
  * (26) V. Balasubramanian, P. Caputa, J. Magan, and Q. Wu, Quantum chaos and the complexity of spread of states, arXiv:2202.06957.
  * (27) T. Ali, A. Bhattacharyya, S. S. Haque, E. H. Kim, N. Moynihan, and J. Murugan, Chaos and Complexity in Quantum Mechanics, Phys. Rev. D 101 (2020), no. 2 026021, [arXiv:1905.13534].
  * (28) A. Bhattacharyya, W. Chemissany, S. Shajidul Haque, and B. Yan, Towards the web of quantum chaos diagnostics, Eur. Phys. J. C 82 (2022), no. 1 87, [arXiv:1909.01894].
  * (29) A. Bhattacharyya, W. Chemissany, S. S. Haque, J. Murugan, and B. Yan, The Multi-faceted Inverted Harmonic Oscillator: Chaos and Complexity, SciPost Phys. Core 4 (2021) 002, [arXiv:2007.01232].
  * (30) A. Bhattacharyya, S. S. Haque, and E. H. Kim, Complexity from the reduced density matrix: a new diagnostic for chaos, JHEP 10 (2021) 028, [arXiv:2011.04705].
  * (31) A. I. Larkin and Y. N. Ovchinnikov, Quasiclassical method in the theory of superconductivity, Journal of Experimental and Theoretical Physics (1969).
  * (32) E. Dyer and G. Gur-Ari, 2D CFT Partition Functions at Late Times, JHEP 08 (2017) 075, [arXiv:1611.04592].
  * (33) J. Liu, Spectral form factors and late time quantum chaos, Phys. Rev. D 98 (2018), no. 8 086026, [arXiv:1806.05316].
  * (34) A. M. Perelomov, Coherent states for arbitrary Lie group, Communications in Mathematical Physics 26 (1972), no. 3 222 – 236.
  * (35) J.-P. Gazeau, Coherent states in quantum physics.  2009\.
  * (36) K. Kawabata, R. Kobayashi, N. Wu, and H. Katsura, Exact zero modes in twisted Kitaev chains, Phys. Rev. B 95 (2017), no. 19 195140, [arXiv:1702.00197].
  * (37) C. Lanczos, An iteration method for the solution of the eigenvalue problem of linear differential and integral operators, Journal of research of the National Bureau of Standards 45 (1950) 255–282.
  * (38) V. Viswanath and G. Müller, The recursion method: application to many-body dynamics, vol. 23.  Springer Science & Business Media, 2008.
  * (39) M. Heyl, Dynamical quantum phase transitions: a review, Rept. Prog. Phys. 81 (2018), no. 5 054001, [arXiv:1709.07461].
  * (40) D. E. Parker, X. Cao, A. Avdoshkin, T. Scaffidi, and E. Altman, A Universal Operator Growth Hypothesis, Phys. Rev. X 9 (2019), no. 4 041017, [arXiv:1812.08657].
  * (41) J. L. F. Barbón, E. Rabinovici, R. Shir, and R. Sinha, On The Evolution Of Operator Complexity Beyond Scrambling, JHEP 10 (2019) 264, [arXiv:1907.05393].
  * (42) J. M. Magán and J. Simón, On operator growth and emergent Poincaré symmetries, JHEP 05 (2020) 071, [arXiv:2002.03865].
  * (43) E. Rabinovici, A. Sánchez-Garrido, R. Shir, and J. Sonner, Operator complexity: a journey to the edge of Krylov space, JHEP 06 (2021) 062, [arXiv:2009.01862].
  * (44) S.-K. Jian, B. Swingle, and Z.-Y. Xian, Complexity growth of operators in the SYK model and in JT gravity, JHEP 03 (2021) 014, [arXiv:2008.12274].
  * (45) C. Yin and A. Lucas, Quantum operator growth bounds for kicked tops and semiclassical spin chains, Phys. Rev. A 103 (2021), no. 4 042414, [arXiv:2010.06592].
  * (46) A. Dymarsky and A. Gorsky, Quantum chaos as delocalization in Krylov space, Phys. Rev. B 102 (2020), no. 8 085137, [arXiv:1912.12227].
  * (47) A. Dymarsky and M. Smolkin, Krylov complexity in conformal field theory, Phys. Rev. D 104 (2021), no. 8 L081702, [arXiv:2104.09514].
  * (48) P. Caputa and S. Datta, Operator growth in 2d CFT, JHEP 12 (2021) 188, [arXiv:2110.10519].
  * (49) M. Carrega, J. Kim, and D. Rosa, Unveiling Operator Growth Using Spin Correlation Functions, Entropy 23 (2021), no. 5 587, [arXiv:2007.03551].
  * (50) J. Kim, J. Murugan, J. Olle, and D. Rosa, Operator delocalization in quantum networks, Phys. Rev. A 105 (2022), no. 1 L010201, [arXiv:2109.05301].
  * (51) A. Kar, L. Lamprou, M. Rozali, and J. Sully, Random matrix theory for complexity growth and black hole interiors, JHEP 01 (2022) 016, [arXiv:2106.02046].
  * (52) E. Rabinovici, A. Sánchez-Garrido, R. Shir, and J. Sonner, Krylov complexity from integrability to chaos, JHEP 07 (2022) 151, [arXiv:2207.07701].
  * (53) F. B. Trigueros and C.-J. Lin, Krylov complexity of many-body localization: Operator localization in Krylov basis, arXiv:2112.04722.
  * (54) E. Rabinovici, A. Sánchez-Garrido, R. Shir, and J. Sonner, Krylov localization and suppression of complexity, JHEP 03 (2022) 211, [arXiv:2112.12128].
  * (55) N. Hörnedal, N. Carabba, A. S. Matsoukas-Roubeas, and A. del Campo, Ultimate Physical Limits to the Growth of Operator Complexity, arXiv:2202.05006.
  * (56) B. Bhattacharjee, X. Cao, P. Nandy, and T. Pathak, Krylov complexity in saddle-dominated scrambling, JHEP 05 (2022) 174, [arXiv:2203.03534].
  * (57) Z.-Y. Fan, Universal relation for operator complexity, Phys. Rev. A 105 (2022) 062210, [arXiv:2202.07220].
  * (58) K. Adhikari, S. Choudhury, and A. Roy, ${\cal K}$rylov ${\cal C}$omplexity in ${\cal Q}$uantum ${\cal F}$ield ${\cal T}$heory, arXiv:2204.02250.
  * (59) W. Mück and Y. Yang, Krylov complexity and orthogonal polynomials, arXiv:2205.12815.
  * (60) D. Patramanis, Probing the entanglement of operator growth, PTEP 2022 (2022), no. 6 063A01, [arXiv:2111.03424].
  * (61) P. Caputa, J. M. Magan, and D. Patramanis, Geometry of Krylov complexity, Phys. Rev. Res. 4 (2022), no. 1 013041, [arXiv:2109.03824].
  * (62) B. Bhattacharjee, S. Sur, and P. Nandy, Probing quantum scars and weak ergodicity-breaking through quantum complexity, arXiv:2208.05503.
  * (63) S. Bandyopadhyay, S. Bhattacharjee, and D. Sen, Driven quantum many-body systems and out-of-equilibrium topology, J. Phys. Condens. Matter 33 (2021), no. 39 393001, [arXiv:2103.02279].
  * (64) A. Bhattacharyya, S. Das, S. Shajidul Haque, and B. Underwood, Cosmological Complexity, Phys. Rev. D 101 (2020), no. 10 106020, [arXiv:2001.08664].
  * (65) A. Bhattacharyya, S. Das, S. S. Haque, and B. Underwood, Rise of cosmological complexity: Saturation of growth and chaos, Phys. Rev. Res. 2 (2020), no. 3 033273, [arXiv:2005.10854].
  * (66) S. S. Haque, C. Jana, and B. Underwood, Saturation of thermal complexity of purification, JHEP 01 (2022) 159, [arXiv:2107.08969].
  * (67) S. S. Haque, C. Jana, and B. Underwood, Operator Complexity for Quantum Scalar Fields and Cosmological Perturbations, arXiv:2110.08356.
  * (68) J. M. Kosterlitz and D. J. Thouless, Ordering, metastability and phase transitions in two-dimensional systems, Journal of Physics C: Solid State Physics 6 (apr, 1973) 1181–1203.
  * (69) T. Senthil, A. Vishwanath, L. Balents, S. Sachdev, and M. P. A. Fisher, Deconfined quantum critical points, Science 303 (mar, 2004) 1490–1494.

##  Appendix A $\boldsymbol{SU(2)}$ coherent states

Following Caputa:2021sib  we present here an explicit identification of the
K-complexity for $SU(2)$ coherent states with an $SU(2)$ operator expectation
value. Similar results may be obtained for $SU(1,1)$ and Glauber coherent
states. Given a Liouvillian

$L=\alpha J_{+}+\alpha^{*}J_{-}+\gamma J_{0},$ (46)

time-evolving the $SU(2)$ lowest weight state we have

$e^{iLt}|j,-j\rangle=\mathcal{N}e^{z(\alpha,\alpha^{*},\gamma)J_{+}}|j,-j\rangle\equiv\frac{1}{\sqrt{(\bar{z}|z)}}|z),$
(47)

by using, for example, Baker-Campbell-Haussdorff formulae. The number
$\mathcal{N}$ is the normalisation. The use of coherent states make the Krylov
basis immediately apparent as

$|O_{n}\rangle=\frac{1}{\sqrt{\langle
j,-j|J_{-}^{n}J_{+}^{n}|j,-j\rangle}}J_{+}^{n}|j,-j\rangle,$ (48)

so that the Krylov complexity operator is

$\hat{K}_{n}=\sum_{n=0}n\frac{J_{+}^{n}|j,-j\rangle\langle
j,-j|J_{-}^{n}}{\langle j,-j|J_{-}^{n}J_{+}^{n}|j,-j\rangle}.$ (49)

For the $SU(2)$ coherent state the Krylov complexity is thus

$\displaystyle C$ $\displaystyle=$
$\displaystyle\frac{1}{(\bar{z}|z)}(\bar{z}|\hat{K}_{n}|z),$ (50)
$\displaystyle=$
$\displaystyle\sum_{n,m,m^{\prime}}\frac{n}{(m!)(m^{\prime}!)}\frac{\langle\psi_{0}|\bar{z}^{m}L_{-}^{m}L_{+}^{n}|\psi_{0}\rangle\langle\psi_{0}|L_{-}^{n}z^{m^{\prime}}L_{+}^{m^{\prime}}|\psi_{0}\rangle}{(\bar{z}|z)\langle\psi_{0}|L_{-}^{n}L_{+}^{n}|\psi_{0}\rangle},$
$\displaystyle=$
$\displaystyle\sum_{n}\frac{n}{(n!)^{2}}\bar{z}^{n}z^{n}\frac{\langle\psi_{0}|L_{-}^{n}L_{+}^{n}|\psi_{0}\rangle}{(\bar{z}|z)},$
$\displaystyle=$ $\displaystyle z\partial_{z}\log(\bar{z}|z),$

which is exactly equal to the normalised expectation value

$C=\langle\bar{z}|J_{0}+j|z\rangle.$ (51)

##  Appendix B Twisted boundary Kitaev chain

In this appendix we would like to present some preliminary results related to
the twisted Kitaev chain. The Hamiltonian and site-dependent chemical
potential of this model are given already in (22). In order to compute the
spread complexity for this model we first need to compute the Krylov basis.
For general parameters this is a non-trivial task. The boundary terms present
in (22) break translational invariance so that a Fourier transform no longer
yields a tensor product structure. As a result of this the eigenstates for the
twisted Kitaev chain can no longer be expressed as $SU(2)$ coherent states.
Obtaining the Krylov basis is thus a more formidable task which we will
circumvent by making a suitable approximation.

The approximation we apply is based on the following considerations. In
Appendix A we showed that the $SU(2)$ K-complexity is, up to a constant, the
expectation value of $J_{0}$ with respect to an $SU(2)$ coherent state. When
performing the inverse Fourier transform on (17) we thus find that the $SU(2)$
K-complexity is the expectation value of the number operator with respect to
this coherent state. Physically the number operator expectation value makes
sense as a good measure of complexity since we need to add fermions to the
vacuum in order to obtain the target states and adding more fermions is
proportionally more costly in terms of circuit operations.

In light of this result we calculate the expectation value of the number
operator as an estimate of the spread complexity. We stress that this is not
the spread complexity for the twisted boundary condition, since we have not
computed the Krylov basis corresponding to the twisted Kitaev chain
Hamiltonian. We expect, however, that the number operator expectation value
will capture similar features. Explicitly, we compute the ground state of the
twisted Kitaev chain and then compute the expectation value

$N=\langle\psi_{gs}|\sum_{i=1}^{L}c_{i}^{\dagger}c_{i}|\psi_{gs}\rangle$ (52)

as our rough estimate of complexity. Since we are working at finite $L$ we
have to, in addition to computing the expectation value, compute the spectrum
numerically and determine whether there is gap closing. As a cross-check of
the accuracy of our numerical results we have compared with the analytical
results of Kawabata:2017zsb  in the cases where exact results for finite $L$
are known.

There are many parameters to adjust in the twisted Kitaev chain. As such, our
numerical studies probably do not cover all interesting features. We do,
however, report the following observation relevant to topological phase
transitions: Over a significant sample size of cases we have found that,
typically, a closing energy gap between the ground state and first excited
state (or equivalently, a degenerate ground state) can be identified by the
expectation value (52). Specifically, these may be identified as
discontinuities in the expectation value in the presence of the twisted
boundary.

As an instructive example, consider the case $J=\Delta=1$ and $a=b=1$. For
these parameters and even $L$ we find, numerically, a simple relation for the
existence of gap closing namely $|\phi_{1}|<|\phi_{2}|$ in the interval
$\phi_{i}\in[-\pi,\pi)$. By this we mean that there exist values for the
chemical potential in the interval $|\mu|<1$ for which the ground state
becomes degenerate. Furthermore, for even $L$ we always observe an even number
of points for which this occurs, see Fig. (7(b)) for $L=8$ as a representative
example of this.

(a)

(b)

Figure 7: The energy gap between the ground state and first excited as a
function of $\mu$ with $\Delta=J=1$, $a=b=1$. We have fixed
$\phi_{1}=-\frac{\pi}{4}$. a) The chain length is $L=8$. For these parameters
we note a pair of values for $\mu$ exhibit a closing energy gap provided
$|\phi_{2}|>|\phi_{1}|$ b) The chain length is $L=9$. In this case we observe
only a single value for $\mu$ for which gap closing occurs.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 8: The expectation value of the number operator as a function of $\mu$
with $\Delta=J=1$, $a=b=1$. We have fixed $\phi_{1}=-\frac{\pi}{4}$. Whenever
the energy gap between the first and excited state closes we observe a
discontinuity in the expectation value. The values of $\mu$ where this happens
are pointed out by a vertical black line which lines up exactly with the
discontinuity.

For odd $L$ we usually find an odd number of values of $\mu$ where this occurs
and it can also occur for some values where $|\phi_{2}|<|\phi_{1}|$. For both
even and odd $L$ cases, see Fig. (8), we observe the advertised discontinuity
whenever the ground state becomes degenerate. The larger the difference
$|\phi_{2}|-|\phi_{1}|$ the bigger the discontinuity.

As a second example consider the case $J=1,\Delta=\frac{1}{2}$ and $a=1,b=2$.
We highlight this case for a two reasons. Firstly, in the case of even $L$ we
no longer find the simple relation $|\phi_{1}|<|\phi_{2}|$ to predict the
existence of gap closing at some value of $\mu$. Secondly, we see the
appearance of a third location where this occurs in Fig. (10) (b). In general
we find an even (odd) number of values of $\mu$ for which gap closing can
occur for even (odd) $L$. Both these observation resonate with the analytical
results of Kawabata:2017zsb  where an even (odd) number of zero modes can be
found for a chain with an even (odd) number of sites and, in general, the
equations determining the existence of zero modes are involved for finite $L$.

The intuition for the observed discontinuity in (52) is rather
straightforward. As we are changing the value of $\mu$ continuously the
spectrum and eigenstates change continuously. Focusing on one of these states,
the expectation value of the number operator w.r.t this state will also change
continuously. However, when we encounter a value of $\mu$ where the ground
state energy gap closes, there is a discontinuous change in the ground state
wave function since the ground state and first excited state wave functions
are exchanged in the spectrum. This accounts for the discontinuities we
observe in the number operator expectation value.

(a)

(b)

Figure 9:  The energy gap between the ground state and first excited as a
function of $\mu$ with $J=1,\Delta=\frac{1}{2}$, $a=1,b=2$. We have fixed
$\phi_{1}=-\frac{\pi}{4}$. a) The chain length is $L=8$. For these parameters
we note a pair of values for $\mu$ exhibit a closing energy gap provided
$|\phi_{2}|>|\phi_{1}|$ b) The chain length is $L=7$. In this case we observe
only a single value for $\mu$ for which gap closing occurs for
$|\phi_{2}|>|\phi_{1}|$. We find three values for which this occurs for
$|\phi_{2}|<|\phi_{1}|$

(a)

(b)

(c)

(d)

(e)

(f)

Figure 10: The expectation value of the number operator as a function of $\mu$
with $J=1,\Delta=\frac{1}{2}$, $a=1,b=2$. We have fixed
$\phi_{1}=-\frac{\pi}{4}$. Whenever the energy gap between the first and
excited state closes we observe a discontinuity develop in the expectation
value. The values of $\mu$ where this happens are pointed out by a vertical
black line which lines up exactly with the discontinuity.

To summarize; in the context of the twisted Kitaev chain the sensitivity of
spread complexity to the topological phase transition is related to the
sensitivity of the number operator expectation value to the ground state
degeneracy. A degenerate ground state signals a discontinuous change in the
expectation value, albeit at finite $L$, due to a discontinuous change in the
ground state wavefunction. The spread complexity, in the large $L$ limit, is
precisely the expectation value of the number operator with respect to some
eigenstate of the Hamiltonian. At finite $L$ this connection is not exact,
since the Krylov basis should be updated in the presence of the twisted
boundaries. Nevertheless, we expect this connection to be a good approximation
to the true spread complexity.