In 2012, Kholodenko and Silagadze published an enigmatic paper¹ in which they showed for n = 1, 2, 3, … that

Their approach involved an intricate use of the concepts from quantum physics, namely, spinorization and the Hopf map, using which they were able to show that the integral is equivalent to the Landau–Zener formula. The Landau–Zener formula expresses the exact probability amplitude of a certain quantum system with a dynamical potential and its standard derivation uses either the contour integration (Landau’s approach) or asymptotics of a certain differential equation (Zener’s approach); see the work of Piquer i Méndez² for a clear exposition of both approaches. Alternatively, see the work of Rojo³ for a matrix exponential solution or that of Glasbrenner and Schleich⁴ for a derivation using Markov property.

In our paper, we define four auxiliary functions of a real parameter t, which turn out to be relevant,

Those functions can be expressed in an exact form using Appell and generalized hypergeometric functions. The exact expressions are, however, not necessary for our derivation of I_n. Note that those functions are not independent. By substituting x → 1/x, we obtain an obvious relation between U and V,

Moreover, the complex-valued functions P, Q can be expressed as linear combinations of real-valued-only functions U and V, respectively. In order to find those relations, we need to study a certain contour integral.

Consider the following contour integral with parameters ⁠, and ⁠:

where we define ⁠. Both  and ln z functions are assumed to have branch cuts at the negative real axis (the so called principal branch), that is, arg  z ∈ (−π, π]. Before proceeding further, let us clarify all the functions appearing in the contour integral. Readers familiar with complex functions may skip the following paragraphs. Writing z = re^iφ with r > 0 and φ ∈ (−π, π], the complex square root function is defined by

so  in the formula above is given by

Similarly, we have for the complex logarithm,

Those definitions ensure both  and ln z are analytic on ⁠.

Note that we recover the ordinary arctan function on reals. To see this, let ⁠, then  and thus

Substituting into the complex definition of arctan,

since 2 arctan x ∈ (−π, π).

To locate the branch cuts of the complex function arctan z, we need to find which points on the z-plane are mapped onto the negative real axis in  as ⁠. Solving

we get the location of the branch cut in  as

where (z₁, z₂) denotes a complex interval from z₁ to z₂ (a line segment between any ⁠).

Let z = iy, y > 0, so z = ye^πi/2 and thus

The situation is different for the arctan function since it is no longer analytic on the whole imaginary axis. Let z = iy, y ∈ (0, 1) (analytic part), then

where  is the usual inverse hyperbolic tangent function. As we have derived, the branch cut of arctan z function is located at the imaginary axis. Let us study how the function behaves in its neighbourhood. Let z = iy + ɛ, y > 1 with ɛ > 0 small positive (so z is taken to lie to the right from the branch cut). By Taylor expansion,

which gives

Let |z| → ∞, expanding ⁠, we get

If  only (that is z → ∞e^iφ with ⁠), we get arg(1 − 2i/z) ∈ (−π, 0) and thus

from which, dividing by 2i,

Let us examine the behavior of arctan z near its singularity z = i in the first quadrant in ⁠. Let z = i − iɛe^iφ, ɛ > 0 and φ ∈ (0, π). Then,

from which

Next, the square brackets part of f(z) has a finite limit as z → i. By L’Hospitals rule,

similarly for  as z → i.

We are now ready to examine integral (4). It is convenient to denote

Let us consider the counterclockwise contour C consisting of integration curves C_k, k = 1, …, 5,

as shown in Fig. 1.

Inside C, function f is analytic, so by the Cauchy integral formula,

Let us now parameterize the integrals on the individual integration curves; we denote  and ⊝ C_k as negatively oriented C_k.

• C₁: z = x,   x ∈ (0, R),   R → ∞,   dz = dx,

  (25)

  As R → ∞, we may express the limit using the auxiliary functions defined earlier,

  (26)

  The integral on the right is trivial, we get

  (27)

• ⁠,

  (28)

  from which J₂ → 0 as R → ∞.

• ⊝ C₃: z = iy + ɛ,  y ∈ (1 + ɛ, R),   dz = idy,

  (29)

  Substituting y → 1/y and letting R → ∞ and ɛ → 0⁺, we get, using auxiliary functions,

  (30)

• ⊝ C₄: z = i − iɛe^iφ, φ ∈ (0, π), ɛ → 0⁺, dz = ɛe^iφdφ,

  (31)

  Since ⁠, J₄ vanishes as ɛ → 0⁺.

• ⊝ C₅: z = iy, y ∈ (1, 1 − ɛ), ɛ → 0⁺, dz = idy,

  (32)

  As ɛ → 0⁺, we get, using auxiliary functions,

  (33)

By the Cauchy integral formula, J₁ + ⋯ + J₅ = 0. Taking the limit as R → ∞ and ɛ → 0⁺, we get for any ⁠,

Comparing the terms with the coefficient A and B, respectively, we get

Solving for P and Q, we get, finally,

Relation (1) can be written recursively by introducing new functions τ_n as

with initial condition τ₀(x) = 1. We have

Let us introduce a new parameter  using which we define a generating function T as

Equation (37) then gets transformed into

which is an integral equation for T(x, t) subject to the boundary condition,

The solution of (40) at x → ∞ completely solves our problem of finding I_n since by expanding T(∞, t) in t when t → 0, we get

Note that rescaling by  with s > 0 and letting s → ∞, we get the stretching relations of T(x, t),

we will see their importance later.

It is convenient to denote

Differentiating the integral Eq. (40) with respect to x, we get

At x = 0, we get ⁠. Second differentiation with respect of x gives

At x = 0, we get ⁠. Yet, another differentiation yields

which gives  at x = 0. One final differentiation results in the following differential equation:

This is an ordinary differential equation (ODE) with parameter t. Its solution at x → ∞ completely solves our problem of finding I_n (expanding T(∞, t) in t). For completeness, the following is the complete list of boundary conditions on T(x, t):

Using the Mathematica  command, we can see numerically how the solution looks like for various parameters t (Fig. 2).

In what follows, we establish a relation between T(∞, t) and T(0, t). Let

Since T(x, t) is continuous, bounded, and having a finite limit T(∞, t), the integral on the right converges for any  with ⁠. Assuming s is a positive real number, we get by substitution ⁠,

It is thus convenient to define ⁠. Imposing the stretching relations (43),

Differentiating e(s, t) with respect to s and integrating by parts with respect to x yields

or in terms of E(s, t),

Writing out the right-hand side using Eq. (45) and changing the order of integration,

The inner integral is trivial and is equal to ⁠. Hence, we get a differential equation for E(s, t),

The solution of this ODE with the appropriate boundary condition  must be

The remaining condition on E(0⁺, t) Eq. (52), therefore, yields as s → 0⁺,

In this section, we establish a relation between T(−∞, t) = 1 and T(0, t). Let

That is, the integration domain is (−∞, 0). Again, the integral converges for any  with ⁠. Assuming s is a positive real number, we get by substitution ⁠,

It is thus convenient to define analogously ⁠. Imposing the stretching relations (43),

Differentiating d(s, t) with respect to s and integrating by parts with respect to x yields

or in terms of D(s, t),

Writing out the right-hand side using Eq. (45) and changing the order of integration,

The inner integral is again trivial and is equal to

Hence, using the definition of α(t) and β(t) Eq. (44), we get a differential equation for D(s, t),

The solution of this ODE with the boundary condition  is

The remaining condition on D(∞, t) Eq. (61), therefore, yields plugging s → ∞,

where U and V are defined by Eq. (2).

The key to our problem is the realization that the solution we found for D(s, t) is valid for all  for which ⁠. The functions α(t) and β(t) are not chosen arbitrarily since they are defined using T(x, t). We may also express them Eq. (44) in terms of d(s, t) as

The value d(±i, t) must be viewed as a limit of d(s, t) for s → ±i (since ⁠). This gives us the final boundary condition. Considering only the case s → i, we get in terms of ⁠,

On the other hand, D(i, t) can be calculated directly from the solution we found for D(s, t). Let s = iσ, σ ∈ (0, 1), then by Eq. (68),

Writing

and integrating out the 1/(1 − y²) term, we get

Note that the integral is now finite as σ → 1. Crucially, as σ approaches 1, since the first part already yields the correct limit ⁠, the term in the square bracket must vanish, otherwise, the limit as a whole would not exist. Hence, we get the following condition:

where P and Q are defined by Eq. (2). Substituting for P and Q from Eq. (36), we get

Comparing the real and imaginary parts, we obtain

which is a linear system of equations for α(t) and β(t), with the solution given as

Substituting α(t) and β(t) into Eq. (69),

from which, by Eq. (58),

The series expansion of T(∞, t) Eq. (42) immediately yields

A similar analysis also shows that