Continuum limit of random matrix products in statistical mechanics of disordered systems
Francis Comets Giambattista Giacomin Rafael L. Greenblatt
math-ph math.MP math.PR 82B44, 60K37, 82B27, 60K35







Continuum limit of statistical mechanics random matrix products]Continuum limit of random matrix products 
 in statistical mechanics of disordered systems


F. Comets, G. Giacomin and R. L. Greenblatt]Francis Comets, Giambattista Giacomin and Rafael L. Greenblatt
[FC,GG]Université Paris Diderot, Sorbonne Paris Cité,  Laboratoire de Probabilités, Statistiques  et Modélisation, UMR 8001,
            F-75205 Paris, France
[RLG]Università Roma 3, Dipartimento di Matematica e Fisica, 
Largo San Leonardo Murialdo 1,
00146, Roma, Italy


We consider a particular weak disorder limit (continuum limit) of matrix products that arise in the analysis of  disordered statistical mechanics systems, with a particular focus on random transfer matrices.
The limit  system is a diffusion model for which the leading Lyapunov exponent can be expressed explicitly in terms of modified Bessel functions, a formula that appears in the physical literature on these disordered systems.  We provide an analysis of the diffusion system as well as of the link with the matrix products. 
We then apply the results to the framework considered by  Derrida and Hilhorst in <cit.>, which 
deals in particular with the strong interaction limit for disordered Ising model in one dimension and that identifies 
a singular behavior of the Lyapunov exponent (of the transfer matrix),
and to the two dimensional Ising model with columnar disorder (McCoy-Wu model). We show 
that the continuum limit sharply captures the  Derrida and Hilhorst singularity. 
Moreover we revisit the analysis by McCoy and Wu <cit.> and remark that it can be interpreted in terms of  the continuum limit approximation. We provide a mathematical  analysis of the continuum approximation of 
the free energy of the McCoy-Wu model, clarifying the prediction (by McCoy and Wu)
that, in  this approximation, the free energy of the two dimensional Ising model with columnar disorder is C^∞ but not analytic at the critical temperature. 


AMS  subject classification (2010 MSC):
82B44, 
60K37,  
82B27, 
60K35 


Keywords: disordered systems, Lyapunov exponents, weak disorder, continuum limit, critical behavior,  two dimensional Ising model, columnar disorder.


[
    [
    December 30, 2023
=====================





§ INTRODUCTION

Products of random matrices can often be interpreted, in a statistical mechanics perspective, as models of disordered systems. 
The leading Lyapunov 
exponent  may then be identified with some physical quantity such as the free energy density  or  persistence length.
We can also take the opposite viewpoint and ask whether a disordered system can be written in terms of, or at least approximated by, a suitable product of random matrices.
 It turns out that there are several examples in which this can be done. Examples include essentially all statistical mechanics systems in which there is a natural one dimensional structure, but it goes also beyond this: the literature is too wide to be properly cited here and we refer to the reviews <cit.>.
 Of particular interest for us are the examples arising from the transfer matrix approach in the statistical mechanics of disordered systems.
 For one  dimensional (let us say, Ising or Potts) models  with finite range interaction   one can write the partition function in terms 
 of a product of matrices <cit.>: 
  if the interactions are only one body and nearest neighbor two body  the transfer matrix of an Ising model is  a two by two matrix, and the size  is larger for Potts and/or longer range models.  But for  two or more dimensional systems
 the size of the  transfer matrix   tends to infinity in the thermodynamic limit and the transfer matrix should be thought more as a transfer operator: this is true also in one dimension if the spin variable can take an infinite number of values <cit.>.
 Nevertheless, also in these cases finite dimensional matrix models can  be helpful (for numerical approximations for example, but also for rigorous bounds, see for example <cit.> and references therein). It is however  remarkable that also the solution of the two dimensional Ising model with nearest neighbor interactions and no external field can ultimately be  expressed in terms of products of  two by two  matrices: this is the essence of several formulations of the celebrated solution of Lars Onsager <cit.>. What is even more remarkable from our viewpoint is that this structure still holds when  special types of disorder are introduced, giving a product of random matrices <cit.>. 
  
  The two by two matrices that arise in the problems we have just mentioned have a particular form: it is

    [ 1  ; Z Z ]  ,

 where  is a real number – say ||≤ 1/2 to keep far from the zero determinant case = ± 1 –  and Z is a positive random variable with log_+(Z) < ∞.  Let us call (informally for the moment)
 _Z() the Lyapunov exponent of a product of IID matrices of the form (<ref>), which appear notably in the following two contexts.

  
 
  * In the one dimensional Ising model with random  external field h=h_j – that is, {h_j}_j=1,2, … is a sequence of independent identically distributed (IID) random variables – and nearest neighbor interaction J, 
 the transfer matrix can be cast in the form (<ref>), with Z=exp(-2h) and =exp(-2J). The free energy density 
 is therefore precisely _Z() 
  and the ↘ 0 limit is the limit of strong ferromagnetic 
 interaction. 
 
  * In a much less straightforward way (detailed in Appendix <ref>), also the free energy of the two dimensional Ising model with a special type   of disordered nearest neighbor interactions (columnar disorder), and no external field, is (essentially) just  ∫_0^1/2_Z(), of course with a proper choice of Z=Z_ that contains  the inverse temperature  of the system.
 The phase diagram of this model (that is, the presence and nature of phase transitions) is determined by the regularity of this expression as a function of ; 
 the most notable prediction for this model, which goes now under the name of McCoy-Wu model, is that the second order transition
 of the two dimensional non-disordered Ising model (for which the second derivative of the free energy diverges
 at criticality like -log| -_c|) becomes of infinite order when the columnar disorder is introduced:  that is, the free energy at the critical point is C^∞ but not analytic. The precise nature of the singularity is characterized in <cit.> by means of a divergent power series for the free energy at _c, where the value of _c depends on the disorder: a summary of the expected effect of disorder on the transition for the two dimensional Ising model is in <cit.>. The McCoy-Wu model has a prominent role in physics because it can be mapped  to the one dimensional quantum spin chain with transversal magnetic field  <cit.>
 and because it
  has played a central role in the development of the real space strong/infinite disorder renormalization group
(see e.g.  <cit.> and <cit.>).
 
 
 Other contexts in which  (<ref>) and  _Z() arise include one dimensional random walk in random environment and  a number of random hopping problems (see <cit.> and references therein), and the key issue for us is that 
 all this vast literature focuses on the → 0 behavior of _Z(), see notably <cit.>.
 From a mathematical viewpoint this limit is of interest because, thanks to <cit.> (see also <cit.>), we know that 
 ↦_Z() is real analytic if ||∈ (0,1) under additional mild hypotheses on  Z (for example: (Z=c)=0 for every c). But the regularity at =0 is not obvious, as well as  if there is a singularity at all. And this is precisely the question addressed in <cit.>. In particular _Z()
  is expected to have a fractional or logarithmic  scaling when → 0 under the  frustration hypothesis that (Z>1) and (Z<1) are both positive: this is the case for example of the ^2 singularity  found in <cit.>, and proven mathematically in <cit.>,  and that we will explain in detail in Section <ref>.  
   
  Here we do not address the study of the Lyapunov exponent of products of matrices of the form  (<ref>). Rather we 
 focus on a continuous time model that arises as a diffusive limit of the matrix product (we call it continuum limit). Roughly, the limit is achieved 
 by considering matrices close to the identity:  is replaced by , with ↘ 0 and we consider Z=Z^ that is very concentrated around one: both [Z^]-1 and 
 var(Z^) are of order  . The dynamics will therefore happen on a timescale 1/ and it will be governed by a two dimensional  stochastic differential system. We then  study the leading Lyapunov exponent () of this limit system: we will actually show that  _Z() ∼() for ↘ 0 (we use ∼ for asymptotic equivalence: the ratio of left-hand and right-hand sides converges to one). This limit has been already considered in several works  and even in greater generality (matrices close to the identity: see e.g.<cit.>), but mathematically rigorous results are lacking (with the exception of <cit.>, whose assumptions however exclude the case we treat): 
 the type of results one finds are expansions of the type
 
    _Z^()   =  c_1()  + c_2() ^2 +…  ,

 where of course c_1()=() and expressions or at least procedures to compute the c_j() are given. To be precise, a full expansion like (<ref>) is not expected to hold in general and, even in the cases in which it holds, e.g. <cit.>, and assuming smooth dependence in  of the coefficients of the matrix,  there is to our knowledge no proof that  
 ↦_Z^()∈ C^∞. We point out, however,  a  very special example in <cit.> that has been worked out explicitly 
and for which the Lyapunov exponent is  analytic except at zero 

where it is nonetheless  C^∞ (note also that <cit.> cannot be applied because for =0
the matrix is the identity matrix).   
 
 We will focus only on (): in other words, the continuum limit we consider captures only   the leading order term in (<ref>). 
 The first remarkable fact is that () has an explicit expression in terms of a ratio of modified Bessel functions: we provide a proof of this fact, which has long been known in the physics literature. To our knowledge, it is found for the first time in  <cit.>, and it then reappears in other works and contexts, see for example <cit.> to which we refer also for a comprehensive review of the literature. It is rather surprising that, while 
 a detailed analysis of the → 0 limit  of () is rather 
 straightforward (the case of ∈ [0, 2) is worked out in  <cit.>), a full analysis appears to be lacking, as well as an emphasis on the rather striking fact that  
 the → 0 behavior   of ()  captures all known and conjectured features of the → 0 behavior   of _Z^ (), i.e. for matrix products (without assuming the  disorder to be small).
  In particular, the ^2 singularity  found in <cit.>  is fully present in the continuum limit expression: mathematical results on this issue for _Z^ ()  have been  recently obtained  <cit.>, but the control of the singular term is an open problem for ||≥ 1. 

Turning to  the  McCoy-Wu model, we come back to the fact that this model appears prominently in the physical literature, in part of course because of its exactly solvable character. 
And  the conventional wisdom in the mathematical community appears to be that 
the  McCoy-Wu claims are exact. And this is correct as far as the free energy formula (in terms of the Lyapunov exponents) is concerned. 
The subsequent analysis is less sound: _c is identified via the equation log Z=0 – the random variable Z depends on the inverse temperature  – and this assertion has some grounds at least at a heuristic level, but then one has to show that the free energy 
∫_0^1/2_Z^() is not analytic at =_c.
And this is (ultimately) done by replacing _Z^() with () and this step 
is very weak on mathematical grounds because the McCoy-Wu claim (which provides the motivation for the whole exercise) is for >0 (possibly very small, but non zero): making this step 
rigorous – possibly by controlling the remainder of the series in (<ref>) – appears to be very challenging, and we do not address this in the present work. 
Once this approximation is done, McCoy and Wu are left with  studying  the regularity in  of ∫_0^1/2().
In spite of being a relatively explicit expression, this is still challenging.  
McCoy and Wu  do this by developing the ratio of Bessel functions in the expression for () for  close to _c
–  the dependence in  is in the index
of the Bessel functions  –  and by identifying the leading-order (in magnitude) terms in an expansion of the Bessel functions as the most singular part. 
We provide a proof that ∫_0^1/2() is C^∞ but not analytic at =_c
and that the asymptotic series at _c is qualitatively the one found by McCoy and Wu (up to a multiplicative factor that they lost
when singling out the most singular term; a similar correction was noted by Luck <cit.> in a related model): technically, this is the most demanding part of our contribution. 
We stress however that what we prove does not yield results  on the  transition for the McCoy-Wu model. The challenging gap pointed out just above remains unclosed.
 But we believe that our contribution helps  understanding the true
 content of the remarkable McCoy-Wu analysis. As a side remark: the computation of McCoy and Wu is done 
 for a very special form of the disorder distribution while they affirm that they expect the result to be true in great generality. 
 We work under very general assumptions on the disorder, thus substantiating this   claim. 



We begin by presenting the diffusion system  and its analysis. This is a 
 stochastic dynamical system that is interesting in its own right and we provide 
 a detailed analysis that goes beyond the strict purpose of what has been explained up to now. 
 In particular we  prove a Central Limit Theorem on the fluctuations
of the Lyapunov exponent  for  the diffusive  limit system, with an explicit formula for the variance and an analysis of the singular behavior.
We then provide a proof that the Markov chain associated to the matrix product described above 
does scale to the diffusion system and that to leading order (in )  the Lyapunov exponent of the Markov chain 
is asymptotically proportional to the Lyapunov exponent of the diffusion system, that is
_Z() ∼(). 
The rest of our work focuses on the regularity/singularity properties of () and of the expressions related to it 
that are of physical relevance. 
 
 

 
  



 §.§ The diffusion model and its leading Lyapunov exponent

We consider the  solution to the stochastic (Itô) differential equations

    X_1(t)  =  X_2(t)  t  ,
    
     X_2(t)  = (  X_1(t)+ (1-) ^2/2  X_2(t) )  t +  X_2(t)  B_t  ,

where B_· is a standard Brownian motion, ≠ 0, ∈  and >0. We consider
deterministic   initial condition (X_1(0), X_2(0))∈^2∖{(0,0)}. The case (X_1(0), X_2(0))=(0,0), as well as =0, are excluded because they are atypical and trivial. The system 
(<ref>) is linear  with a multiplicative noise so, given the initial condition, there exists a unique strong solution. Our focus is on the Lyapunov exponent ()=_,  () that we introduce via our first statement in which we use the Euclidean norm ‖·‖ in ^2 just for definiteness. Before stating it we need to recall
 one of the  definitions of the modified Bessel function of 2^nd kind of index ∈ and argument x>0
10.25

    K_ (x)  := ∫_0^∞exp(-x cosh(t)) cosh( t)   t   = 1/2∫_0^∞1/y^1+ exp(- x/2(y + 1/y))
     y  .

 We note from now that
K_ (x)=K_-(x).






For every ≠ 0 and every 
(X_1(0), X_2(0))∈∖{(0,0)} the limit 

    lim_t →∞1/tlog‖ (X_1(t), X_2(t))‖  =:  _,  ()  ,

exists and  does not depend on (X_1(0), X_2(0)).  Moreover


  * the limit is unchanged if we replace ‖ (X_1(t), X_2(t))‖ with | X_j |, j=1,2 as well as if 
we remove the expectation (in this case the convergence is almost sure);

  * if >0 [resp. <0], then (X_1(t))= (X_2(t)) [resp. (X_1(t))≠(X_2(t)))]
for all t≥τ:= inf{t≥ 0:  (X_1(t))= (X_2(t))} [resp. τ:= inf{t≥ 0:  (X_1(t))≠(X_2(t))}]
and [τ ]< ∞. Moreover 
_,  ()=_,  (-).

  * For >0 and every ∈ we have 

    _, ()  = ^2 /4(x K_-1(x)/ K_(x)) ,      with   x:=4 /^2 .








We draw the attention of the reader on the identification of x with 4/  ^2. This shortcut notation is kept in all statements and proofs.






Some of the results in Theorem <ref> can be understood on the basis of a symmetry enjoyed 
by our system: If (X_1(·),X_2(·)) solves (<ref>), then 
(-X_1(·),X_2(·)) solves (<ref>) with  replaced by -. 
We can therefore restrict our analysis to the case >0 and we will show that 
only the (interior of the) quadrants in which both coordinates have the same sign – first and third quadrant – are recurrent for the dynamics:
all the rest is transient. By linearity we can then restrict to the first quadrant.  



As already mentioned in the introduction, most of the content of Theorem <ref> is known in the physical literature
and (<ref>) appears in number of contexts. Besides the 
 pioneering work <cit.> that we have already mentioned,    (<ref>)  appears for example also in <cit.>, which deals with disordered quantum Ising chains with transverse magnetic field: this is not surprising, because 
this quantum model is mapped exactly (by a Suzuki-Trotter path integral) into a suitable limit of
the McCoy-Wu model  
(this is also  exploited   <cit.>). It also appears in the analysis of one dimensional random Schrödinger equation and in   the analysis of a diffusion of a particle in a random force field, see e.g. <cit.>, and <cit.> which is possibly the first work addressing precisely what we refer to as the continuum limit. Mathematical works have also been done in this context 
and  (<ref>) appears in  a study of the quenched large deviations of diffusions in a random environment <cit.>:
By Kotani's formula (unpublished, 1988)
the Laplace transform of hitting times for the diffusion can be expressed in terms of a Riccati equation – the K_(·) Bessel function appears as solution of this equation – that is equivalent  to (<ref>) below. 

Finally, the first  item  of Theorem <ref> is a  classical result at the random matrix level, and the second item is an elementary observation. In the continuum set-up, the first two items follow by applying standard tools of stochastic analysis (the proofs turn out to be rather concise and we give full details). The third item is a computation: it is not novel, but it is very short and we provide it for completeness.  



We also have a rather explicit representation for the fluctuations:


 
      
The family of random variables

    {1/√(t)(  log‖ (X_1(t), X_2(t))‖ - t _,  () ) }_t∈ [0, ∞)

converges in law for t →∞
to a centered Gaussian variable with variance v_,  ()   ∈ (0, ∞),

    v_,  () =
     2/^2 K_(x)∫_0^∞1/y^1-e^x/2(y+1/y)( ∫_0^y  z-_, ()/z^1+ e^-x/2(z+1/z) z )^2  y
     .





In physics   the  behavior of fluctuations for matrix products has beed repeatedly addressed, see <cit.> and references therein; the same is true for the mathematical literature <cit.>. Proposition <ref> is about the fluctuations for the continuum limit: even taking into account the results  in  <cit.>, we think that
understanding how random matrix products fluctuations and continuum limit fluctuations are related is an open issue.
 









 §.§ Small  asymptotic expansions


Thanks to Theorem <ref> item (3), studying the small  behavior of _, () is 
 just a book-keeping exercise that exploits  
the asymptotic behavior of K_(·). 
For  ∈ [0, 2) this result can be found in   <cit.>, see also
<cit.> for =0: 
we provide the general result and 
we  will explain the relevance of this exercise in Section <ref>. 
Throughout this work Γ denotes the  Gamma function, see  5.2 for definitions
and properties.





Recall that x=4 /  ^2.
For ∈ (0,∞) ∖ we have for ↘ 0

    4/^2_, ()   =  
    c_1 ()x^2+ …  + c_⌊⌋ ()x^2⌊⌋ +2
    Γ(1-)/Γ()( x/2)^2 +O(x^min(2⌈⌉, 4 )) ,

where c_j(·) is a rational function (for explicit expressions, see (<ref>)-(<ref>)).

For ∈{1,2, …} we have

    4/^2_, ()   =  
    c_1() x^2+ … + c_ -1()x^2( -1)+ (-1)^2^2-2/((  -1)!)^2
    x^2log x+
    O( x^2)  ,

where  c_j( · ) is the same rational function as in the non integer case. 

For =0 we have 

    _, ()  =  ^2/4 log (1/x) +O ((log 1/x )^-2)  ,

and
the result for <0 is directly recovered from (<ref>)-(<ref>) by using the identity

    4/^2_, () <0= 2 ||+ 4/^2_, ||()  .






The identity (<ref>)
is a simple consequence of the Bessel identity

    x K_1+(x) =  2  K_(x)+ xK_ -1+(x)  ,

that follows from (<ref>) by integration by parts, together with the identity K_α(x) = K_-α(x).




We have chosen to give these expansions up to the leading singular term: one can of course be much more precise. 




Keeping only the leading term, Proposition <ref> implies

    _σ,α() ↘ 0∼(^2/4)
    
    c_1() x^2     if  >1  ,
    
    2Γ(1-)/Γ()(x/2)^2    if ∈(0,1)  ,
    
    
    1/log (1/x)       if = 0  , 
    
    2 ||    if ∈ (-∞,0) ;

recall again that x = 4ε/σ^2.
It is certainly worth observing for the benefit of those readers who are less at home with special function calculations that  
(<ref>) can be derived by elementary asymptotic methods from 
(<ref>). This is of course also the case for the full Proposition <ref>, but the exercise 
becomes particularly involved and resorting to the special functions literature is certainly wise, or  even necessary.






	(<ref>) directly entails lim_↘ 0_, ()=^2/2||_ <0,  and ↦^2/2||_ <0 is singular at the origin, while 
↦_, () is real analytic for ≠ 0 (and it is meromorphic in the whole : see beginning of Section <ref>).  It is possibly worth
observing that we have not defined _, () for =0 because of the pathological nature of this case, but 
the limit in (<ref>) exists also for =0 and  _, (0)=
^2/2||_ <0 (in agreement with lim_↘ 0_, ()), but only if  X(0) ≠ 0; otherwise the Lyapunov exponent is -∞. 
Moreover
the (Laplace) asymptotic behavior of the two components for =0 in general does not coincide with the 
Lyapunov exponent. 




Of course one could wonder about the behavior as ↘ 0 of the variance v_, () in
 Proposition <ref>.  





Still with the notation x=4 /  ^2, we have that for every  there exists C()>0 
(see (<ref>) for an explicit expression) such that

    v_, () ↘ 0∼ C() ^2/2×
     1     if ≤ 0  ,
    
     x^2    if ∈ (0,2)  ,
    
      x^4log (1/x)     if =2  , 
    
     x^4    if >2  .







 

 §.§ From matrix product to the diffusion model

 We are now going provide rigorous results about  how  (<ref>) emerges  as limit of matrix products. 
Consider, for  >0 and given  an IID sequence {_n}_n=1,2, …  of standard Gaussian variables,  the discrete time stochastic 
process {(X_1^(n),X_2^(n))}_n=0,1, … defined recursively from the deterministic initial condition (X_1^(0), X_2^(0))=
(X_1(0), X_2(0)) by

    X_1^(n+1)  =  X_1^(n) +  X_2^(n)   ,
    
    
    X_2^(n+1)  = 
    e^√() _n+1 -  ^2/2( X_2^ (n) +  X_1^(n) )  .

 Defining 

    Z^(n+1)= e^√() _n+1 -  ^2/2,

 we can write (<ref>) as 

    X^ (n+1)= X^(n) + A^ (n+1)  X^  (n),

where

    X^ = [ X_1^; X_2^ ] ,
         A^(n) = 
    [       0        ;   Z^(n) Z^(n)-1 ] 
    .

In different terms: X^(n) results form the product of n independent matrices of the form
I+A^ and  for =1 we have that the matrix I+A^1 coincides 
with
(<ref>) when Z=Z^1, that is when Z is log-normal. The restriction to log-normal is just for ease
of exposition: we are going to prove a result (Theorem <ref>) for much more general distributions. 



Note that the determinant of I+A^ is Z^ (1- ^2 ^2) and we want to exclude the degenerate case:  since
we are going to give a result for  ↘ 0, we can assume that this requirement is automatically satisfied. 
The rate of growth of X^(n) is defined by the Lyapunov exponent

    _Z^ () = lim_n →∞1/nlogX^(n) ,

which exists a.s. and is deterministic, see e.g. <cit.>.





 

For ↘ 0, the random process
 
    {(X_1^(⌊ t/⌋),X_2^(⌊ t/⌋))
    }_t∈[ 0, ∞)  ,

converges in law to the diffusion (X_1(·), X_2(·)) on the Skorokhod space 𝒟([0, ∞), (0, ∞)^2).
Moreover,

    lim_↘ 0_Z^ ()/ = _,  ()  .

 









 §.§ Continuum limits and the Derrida-Hilhorst singularity



In <cit.>, B. Derrida and H. J. Hilhorst study the ↘ 0 limit of the Lyapunov exponent
_Z ()
of product of IID matrices of the form (<ref>), under the hypothesis that 
[Z]>1 and [log Z]<0. Since ↦[Z^] is convex, by the hypotheses on Z there exists a unique ≠ 0 such that Z^ =1, and
one readily realizes  that ∈ (0,1). 
It is claimed in <cit.> that

    _Z () ↘ 0∼ C ^2  ,

with a semi-explicit expression for C=C_Z>0, that depends on the law of Z. Such a result directly implies 
a corresponding result for the case [log Z]>0 and [Z^-1]>1: note that in this case [Z^]=1, ≠ 0,
is again uniquely solved and ∈ (-1,0). So, by writing 

    [ 1  ; Z Z ]  =  Z  [ Z^-1 Z^-1;         1 ] ,

we see that 

    _Z ()   = log Z + _1/Z ()(<ref>)=  log Z + C_1/Z^-2 + o(  ^-2) ,

and we recall that ∈ (-1,0) now.

Moreover one can find in <cit.> an argument telling us that for >1,
∉, one 
expects

    _Z () =  c_1 ^2+ … +… + c_⌊⌋^2⌊⌋ +
    C ^2 + o(^2 )  ,

for  real constants c_j and C that are in principle computable. For the case =0, 
 i.e. for the case log Z=0 in which the only solution to Z^=1 is =0, one finds 
the prediction 

    _Z () ↘ 0∼C/log (1 /)  ,

with C>0, in more than one reference. We mention here
 <cit.> in which (<ref>)
 is found for one dimensional Ising model with random field for a very specific choice of  the disorder (the interaction 
  J of <cit.> corresponds to log (1/)). In the localization context (<ref>) has been found for example in
<cit.>.

From our perspective, the significance of all this is that:




  * the continuum limit results of Proposition <ref> fully match with the expected behaviors (to all orders!),  
	(<ref>),  (<ref>) and (<ref>), for the random matrix product. We consider this to be rather striking, and it 
 highlights the richness of  the continuum limit;
 
  * we are going to review the mathematical results available about  (<ref>),  (<ref>) and (<ref>), but we want to point out that even at the level of physical predictions some results are more sound than others. Notably, it appears to be rather challenging to capture the ^2 singularity for || >1 and the level of sharpness of the  ||∈ (0,1) prediction (<ref>), even leaving aside mathematical rigor, does not appear to be 
	 easy to achieve. In this sense, the continuum limit goes beyond what has been established so far for the discrete case.  
  


From a mathematical standpoint
a proof of (<ref>), and (<ref>), (i.e., (<ref>) with ||∈ (0,1)) has been achieved only recently and 
under the assumption that Z has a C^1 density and that the support of Z is bounded and bounded away from zero <cit.>. It is well known, see e.g. <cit.>, that the problem of computing the Lyapunov exponents boils down to
identifying the invariant probability of a Markov chain associated to the matrix product. The arguments in <cit.>
aim at constructing a probability that for  small is expected to be close to the invariant probability.
In <cit.> this construction is put on rigorous grounds and, above all, it is shown that 
this probability, although not invariant, is sufficiently close to the invariant one to make it possible to control the Lyapunov exponent with the desired precision. A result about (<ref>), i.e. for ||≥ 1, has been achieved recently <cit.>, but the the expansion is fully controlled only up to (and excluding) the singular term C ^2: for the moment results about this term remain very weak.





 §.§ On the two-dimensional Ising model with columnar disorder (McCoy-Wu model)
 

It is possibly somewhat unexpected, but also computing the free energy of 
the two dimensional Ising model with columnar disorder
(McCoy-Wu model <cit.>) boils down to analyzing 
the Lyapunov exponent _Z. 
The McCoy-Wu prediction is remarkable and folklore says that their model is the only non trivial  exactly solvable disordered statistical mechanics model:  we dedicate  Appendix <ref>
to introducing in detail the model, keeping close to the McCoy-Wu notations. But the key point from the result viewpoint is that 
L. Onsager celebrated solution of the non disordered
case establishes that the free energy, as function of the temperature, has a (logarithmic) divergence
in the second derivative 
at the critical  temperature. B. M. McCoy and T. T. Wu predict that if a small amount of columnar disorder
(i.e. one dimensional: vertical bounds couplings are random and they are repeated – i.e. no new randomness is introduced – on each line) is introduced  
the transition persists but disorder is relevant (in the sense of the Harris criterion, that is the disorder changes the critical behavior, see e.g.
<cit.>) and the transition becomes C^∞. A precise form of the singularity is also given. 

As  explained in Appendix <ref>, McCoy and Wu extend Onsager's approach to the columnar disorder case and the free energy can be written, up to additive analytic terms, in terms of an
integral in the  variable of the Lyapunov exponent _Z(), with Z that has an explicit expression in terms of the parameters of the Ising model. The analysis by McCoy and Wu 
of this expression is performed in two steps:




  * They claim that in the limit of very narrow disorder 
the relevant – i.e. singular – contribution to the free energy can be written as 

    : ↦∫_(0,η)  x
     K _-1( x)/ K_ ( x) x  ,

with η>0 arbitrary (the singular part comes from the small x behavior of the integrand). 


The integrand is just 4_1,  (x/4), and so it is clear from the estimates in Proposition <ref> that the integral is well defined for all real .

  * They argue, by  approximating the integrand by another expression for which 
the exact integration can be performed, that (<ref>) is C^∞ but not analytic at =0.



The approximation in the first step, see Appendix <ref>, turns out to be precisely the diffusive limit we deal with:
this was possibly expected by comparing (<ref>) and (<ref>).
What we do with the next result is 
providing     a rigorous analysis of the second step, that is the analysis of (<ref>).





 is real analytic in (-1,1)∖{0}. Moreover it
 is C^∞ but not analytic in 0.  The radius of convergence of its  Taylor series at the origin   ∑_n=0^∞ c_n ^n  is zero: in fact  c_1=4η,  c_2n+1=0 for every n∈ and  the even coefficients satisfy 
 
    c_2n n →∞∼   4 e^-γ (-1)^n+1(2n-1)!/π^2n  ,

with γ  the Euler-Mascheroni constant.




We are going to prove more. Namely that (<ref>) defines an analytic  function for every ∈
with 0<|| <1. We believe that the restriction to  0<|| <1 
can be removed to get simply to || >0. However this involves a certain number of complications 
connected to the fact that, with our approach, an ad hoc analysis has to be developed for ∈. 
Since the focus is on =0, we have made the choice not to develop this issue.  



A (very) substantial gap remains between where our results lead and the proof of the McCoy-Wu claim that the
transition is C^∞, even without the precise claim on the nature of the singularity. What we perform, and 
what McCoy and Wu do, is capturing the behavior the free energy near criticality when the disorder is vanishing – this is  reminiscent 
of intermediate disorder limits <cit.> in which, like for us, one enters the framework of  integrable models  – while the true issue is the behavior for (possibly) weak,  but non vanishing, disorder. 






§ ON THE LYAPUNOV EXPONENT: THE PROOF OF THEOREM <REF>









We use the short-cut notation := ^2(1-)/ 2 ∈.
We recall that we can assume >0 and let us start by showing that
the process does not hit (0,0). Recall that (X_1(0),X_2(0))≠ (0,0) and set τ_(0,0):=
inf{ t>0:   (X_1(t),X_2(t))= (0,0)}.
For this let us consider R(t):= √(X_1^2(t)+X_2^2(t)).
By Itô's formula:

    R(t)    = X_1/R X_1+ X_2/R X_2 + 1/2X_1^2/R^3⟨  X_2,  X_2 ⟩
       = (2 X_1 X_2/R + X_2^2/R +  ^2/2X_1^2X_2^2/R^3) t
    +  X_2^2/R B_t
    
       = 
    R (
    2 Y/1+Y^2 + Y^2/1+Y^2 +
    ^2/2Y^2/(1+Y^2)^2)  t + R ( Y^2/1+Y^2)  B_t
    
       =:  R  D t+ R  Q B_t
     
     ,

where Y:=X_2/X_1∈ [-∞, ∞] and D=D(t) and Q=Q(t) are uniformly bounded continuous stochastic processes (‖ D‖_∞≤ 2  + || +  ^2/2 and Q∈ [0, ]),
defined up to τ_(0,0).
Since, again by Itô's formula, we have 

    log R(t)  = (D(t) -1/2 Q^2(t) ) t + Q(t)  B_t  ,

we  see that R(t)/R(0) is bounded away from zero on every compact time interval. 
This readily yields a contradiction if (τ_(0,0)< ∞)>0. Hence  (τ_(0,0)< ∞)=0
and we have proven that the process does not hit the origin.



Now we are going to show that if the initial condition is in the (interior of the) second or fourth quadrant (in the counterclockwise sense), it hits the boundary of these quadrants in an a.s. finite time (in fact, this random time has finite expectation) and enters either the first or third quadrant.
And we show also that once the process is in the first (or third) quadrant, it stays there forever.  

Without loss of generality let us assume that X_1(0)<0 and X_2(0)>0 (second quadrant). 
For the analysis it  is helpful to consider 
Y(t)(<0) up to t≤τ_-∞:=inf{t≥ 0:   Y(t)=-∞}, which coincides 
a.s. with inf{t≥ 0:   X_1(t)=0}, and up to t ≤τ_0:=inf{t≥ 0:   Y(t)=0}. 
By Itô formula

    Y  = ((1- Y^2 )  +   Y )  t +  Y  B_t  ,


and with the specific initial initial conditions we are using is somewhat helpful to 
work with the positive process Ỹ=-Y:
 
    Ỹ  = (( Ỹ^2 -1 )  +  Ỹ)  t + Ỹ B_t  ,

which is in (0, ∞) as long as the two dimensional process does not leave the interior of the quadrant.
We use the stopping times τ̃_0 and τ̃_∞ with the obvious meaning. 
We are going to apply the Feller test for explosion  to show that  τ̃:=min( τ̃_0, τ̃_∞) is in L^1 so 

    ( τ̃< ∞))  =  1  ,

which means that, almost surely, the process hits the axes. 
And if τ̃_∞< τ̃_0, that is if (X_1(τ̃),X_2(τ̃))=(0,x_2), x_2>0, we readily see
from (<ref>) that X_1(τ̃+t)>0, at least for t>0 small.
If instead τ̃_0< τ̃_∞, then (X_1(τ̃),X_2(τ̃))=(x_1,0), x_1<0, 
and again from (<ref>) one sees that X_2(τ̃+t)<0 for t>0 small: since the equation solved by X_2 is stochastic, the argument is slightly more delicate than for the previous case and we give some details. 
By the Strong Markov property it suffices to consider (X_1(0),X_2(0))=(x_1,0), x_1<0, and 
X_2(t)= ∫_0^t ( X_1(s) +c X_2(s))  s +M(t), with the constant c and the centered Martingale M easily read out of  (<ref>). Note that M is a time changed Brownian motion. By continuity of (X_1(·), X_2(·)) we readily see that X_2(t)-M(t) ≤ -| x_1| t/2 for t small. It is therefore  clearly impossible that
inf{t>0: X_2(t)≠ 0} is positive, because this would imply X_2(t)<-| x_2 | /2 for small t. Therefore t↦∫_0^t X_2(s)^2 s is increasing  at least for t small, which implies that  the time change is non degenerate at least for small times. Hence M(t) becomes negative for arbitrarily small values of t. Therefore   X_2(t) < -t| x_1|/2 – in particular, it is negative – for arbitrarily small values of t.
An application of the Feller test, this time applied to Y and not to Ỹ, actually shows that
if Y is in  (0, ∞), then it will stay so for all times, that is the interior the first and third quadrants are stable sets
for the dynamics.

Let us detail the application of the Feller test. Let Z is a one dimensional diffusion with
Z(0)∈ (0, ∞) and

    Z (t)  =  b(Z(t))  t + q(Z(t)) B_t  ,
 
b(·) and q(·)(>0) differentiable functions. We set τ:= inf{ t>0:   Z(t)=0 or Z(t)=∞} and

    s(z)  := ∫_1 ^z exp( -2 ∫_1^y b(r)/q^2(r) r) y
      and v(z)  := ∫_1^z s'(y)( ∫_1^y 2/s'(r)q^2(r) r )  y  .

By monotonicity the limits of v(z) for z↘ 0 and z ↗∞ exist in [0, ∞]
and they will be simply denoted by v(0) and v(∞). 
If both v(0)<∞ and v(∞)<∞ then 
[τ]< ∞ <cit.>. On the other hand, if v(0)=v(∞)= ∞ then 
(τ= ∞)=1 <cit.>  .  























Let us start with the Ỹ case (cf. (<ref>)): we have 

    s'(z) = C/z^1-exp(-2/ ^2(z+1/z) )
     ,

so for z≥ 1

    s'(z) ≍
    z^-1+exp(-2/ ^2z )   ,

with the notation f(z) ≍ g(z) if 
f(z)/ g(z) ∈ [a, 1/a] on the prescribed interval for some  a ∈ (0,1),

and 

    v'(z) ≍1/z^2  ,

so v(∞)< ∞. In a very similar way, for z ∈ (0, 1]

    s'(z) ≍
    z^-1+exp(-2/ ^21/z)   ,

and the estimate of v(0) is identical to  the one for v(∞)  via a (double) change of variable z ↦ 1/z. Hence 
v(0)< ∞ and the diffusion Ỹ hits 0 or ∞ at a random time  which has finite expectation.





For the case of Y we turn to (<ref>) and the difference is that
the factor (z+1/z) in the exponent in (<ref>) changes sign. Once again, we can replace
(z+1/z) by z for z≥ 1, and by 1/z for z ≤ 1. This implies that the integral with respect to 
r in the expression for v(z) in (<ref>) stays bounded and bounded away from zero
both for y↗∞ and for y ↘ 0. The integral with respect to y therefore diverges 
both for z ↗∞ and z ↘ 0 (once again, the two computations are identical, up to change of variables). Therefore, almost surely, Y hits neither 0 nor ∞.



We are now going to show that the diffusion Y has a unique invariant probability, that we will make explicit,
on (0, ∞).  This corresponds to the two (extremal) invariant probabilities
for the normalized process (X_1, X_2)/ √(X_1^2+X_2^2), supported on the intersection of the unit circle with the 
first (or third) quadrant. For this it is practical 
to observe that  the generator of the evolution (<ref>) acts on C^2 functions f: (0, ∞)→ as

    L_ f(y)  = (  (1-y^2) +   y) f'(y) +^2/2 y^2 f”(y)  = ^2/2 p_(y)( y^2 p_(y) f'(y) )' ,

where p_(·) is the probability density 

    p_(y)
      = C_/y^1+ exp(- 2 /^2(y + 1/y)
    )     with   C_^-1  =  2 K_( 4 /^2) ,

and K_(·) is defined in (<ref>).
This already makes evident the reversible nature of the diffusion Y and, in particular,  
(<ref>) is an invariant probability. The transformation S(t):= log Y(t) makes things even more straightforward:
S is a diffusion on  with constant diffusion coefficient and a strongly confining potential: 

    S  =  -U'(S)  t +  B_t    with  U(s)  := ( exp(-s) +exp(s) - ( - ^2/2) s)  .

An invariant probability of this diffusion is p̃_(s)∝exp(-2 U(s)/^2) and the generator has the familiar
symmetric form L̃_ g= (^2/2) ( p̃_ g')'/p̃_, for g ∈ C^2(, ), see e.g <cit.>.

Uniqueness of this invariant measure as well as ergodic  properties can be established in a variety of ways:    <cit.>
gives a Pointwise Ergodic Theorem  that one can directly apply to (<ref>), and of course it implies uniqueness. Alternatively one
can put (<ref>) or (<ref>) in natural scale via a time change and a scale function, see  <cit.>, and apply 
the Ergodic Theorem <cit.>.  Ergodic properties of S are also given in   <cit.>. 


Therefore for every choice of Y(0)∈ (0, ∞), almost surely  and in L^1 we have that

    lim_t →∞1/tlog X_1(t)     = lim_t →∞1/t∫_0^t Y(s) ds 
       = ∫_0^∞ y   p_ (y)  y 
       =  K_-1(4 /^2)/K_(4 /^2)  ,

where in the first step we have used the first identity in

    X_1(t)     =  X_1(0)exp( 
    ∫_0^t  Y(s)  s)  ,
    
    
    X_2(t)     =  X_2(0)exp( ∫_0^t 1/Y(s) s  -  ^2/2 t +  B_t)  ,

which is directly derived from (<ref>) and holds for all t>0 if both X_1(0) and X_2(0) are positive
(or both are negative: in general,  the formula holds up to the hitting time of the boundary of the quadrant in which (X_1(0), X_2(0)) lies).
The second step in (<ref>) is the application of the Pointwise Ergodic Theorem and the last one is an explicit computation.
In the same way, by using the second identity in (<ref>)
we get to (with x=4 /^2)

    lim_t →∞1/tlog X_2(t)     = lim_t →∞1/t( ∫_0^t 1/Y(s) ds) -  ^2/2
       = ∫_0^∞1/y  p_ (y)  y -  ^2/2
       = ^2/4( x K_1+ (x)/K_ (x) -2 )
     (<ref>)=  
      ^2/4x K_1- (x)/K_ (x)  =  K_-1(4 /^2)/K_(4 /^2)  ,

which coincides with what we found in (<ref>). This shows that both components have the same exponential growth rate, hence also the norm of (X_1(t), X_2(t)), and (<ref>) is proven.
If instead of starting from the first quadrant, we were starting from the second quadrant, the result is unchanged because the second quadrant is abandoned after a random time that is in L^1. This completes the proof of Theorem <ref>.




§ FLUCTUATIONS OF THE LYAPUNOV EXPONENT: PROOFS



a

Proof of Proposition <ref>.
Recall that the ratio Y(t)=X_2(t)/X_1(t) converges.  
Hence it is sufficient to prove the convergence of  (<ref>)
with X_1(t) instead of ‖ (X_1(t), X_2(t))‖ in the logarithm, i.e., to prove convergence in law of

    {1/√(t)( ∫_0^t f(Y(s)) ds  ) 
     }_t ∈ (0, ∞)  ,

where the function f(y)= y- _,  () is centered for the invariant density p_(·). We start by solving the Poisson equation,
L_ g = f.
 <cit.> applies for S=log Y, see (<ref>), and shows that the Poisson equation has g(y)= ∫_0^∞_y[ f(Y(s))] ds as unique solution.
We need here an explicit form, and we solve the linear equation

    ^2/2 y^2 h' + (  (1-y^2) +   y) h   = f

for h=g' by the method of variation of constants. The homogeneous equation – when the  right-hand side of (<ref>) is equal to 0 – admits
h_0(y) = (y^2 p_(y))^-1 as a solution. Looking now for solutions of the form h(y)=k(y) h_0(y) for  (<ref>) itself,
we find that 

    k(y) = 2/^2∫_0^y  ( z-_,  () 
     ) 
     p_(z) dz + C ,

and we choose C=0 (this is the only choice that yields the required integrability properties in what follows).
Finally,

    g'(y)  = (y^2 p_(y))^-12/^2∫_0^y  ( z-_,  () 
     ) 
     p_(z)  z ,

and the value of g(1) does not matter for our purpose. 
Now, we can  follow a standard proof of Central Limit Theorem  for reversible diffusions, e.g. <cit.>. 
Note in fact that g is smooth 
and that for y →∞ (recall the notation used in (<ref>))

    g'(y)  ≍  y^ -1e^/2 y(∫_y^∞ z^- e^- /2 z z) 
     ≍ 1/y  ,

where in the first step we use that ∫_0^y  ( z-_,  () 
) 
 p_(z)  z= ∫_y^∞  (_,  () - z
) 
 p_(z)  z and that _,  () is just a constant. For y↘ 0 instead 
 
    g'(y)  ≍  y^ -1e^/2 y( 
    ∫_0^y z^-1- e^-/2 z z )    =  
     y^ -1e^/2 y( 
    ∫_1/y^∞ z^-1+ e^- z/2  z )   ≍  1  .

Therefore sup_y | g'(y) |< ∞ and by  Itô's formula we obtain that

    M_t  := ∫_0^t Y_s g'(Y_s)  B_s     =   g(Y_t)-g(Y_0) - ∫_0^t L_ g (Y_s)  s 
        =  g(Y_t)-g(Y_0) - ∫_0^t f (Y_s)  s    ,

is a martingale with bracket 

    ⟨ M ⟩_t  = ^2 ∫_0^t  Y_s^2 g'(Y_s)^2  s  .

By the ergodic theorem, as t →∞, almost surely

    1/t⟨ M ⟩_t
       ⟶^2  ∫_0^∞ y^2 g'(y)^2 p_(y)  y 
        = 4/^2∫_0^∞1/y^2 p_( ∫_0^y ( z-_,  () 
     ) p_(z)  z )^2  y   
     =  v_,  ()  .

This  deterministic limit is  finite  in view of the (exponential) decay of ∫_0^y ( z-_,  () 
 ) p_(z)  z as y → 0 and y →∞. Then, 
the central limit theorem for martingales applies, and 
 t^-1/2M_t converges in law to a centered Gaussian with variance given by 
v_,  ().  Now, the first two terms in the last line of (<ref>) are bounded in probability, so 
-t^-1/2∫_0^t f(Y_s)ds converges to the same limit as t^-1/2M_t, and (<ref>) is proved. 
Therefore the proof of Proposition <ref> is complete.






Proof of Proposition <ref>. 
 In view of  (<ref>) and of the fact that we know the asymptotic behavior of K_(x)∼ x^-||Γ (||)/2^1-|| for ≠ 0 and K_0(x) ∼log (1/x), what we have to estimate is 
 
    ∫_0^∞1/y^1-e^x/2(y+1/y)( ∫_0^y  z-_, /z^1+ e^-x/2(z+1/z) z )^2  y
      =
    
    ∫_1^∞1/y^1-e^x/2(y+1/y)( ∫_y^∞ z-_, /z^1+ e^-x/2(z+1/z) z )^2  y +                  
    ∫_0^1 1/y^1-e^x/2(y+1/y)( ∫_0^y  z-_, /z^1+ e^-x/2(z+1/z) z )^2  y   =:  T_1(x)+T_2(x)  .






In view of Proposition <ref> we know that   _, =O(^min(2, 2)) for >0, except for =1
for which there is a  logarithmic correction. Therefore _,  =o(  ) if  >1/2 and, since z≥ 1, in dealing with  T_1(x) we can safely neglect the term containing _, for >1/2.  On the other hand, 
in dealing with  T_2(x) we can safely neglect the term not containing _, for <1/2. In fact 
_, is much greater than , hence of z (z ≤ 1 for T_2), for <1/2. 



To make the expressions more compact and readable we choose 

    ^2  =  2  ;

the general case is easily recovered by a scaling argument.




We start with the analysis of T_1. By a change of variable we have:

    T_1(x)   = ( x/2)^∫_x/2^∞ y^ -1e^y+ (x/2)^2/y( ∫_y^∞( z^-  -_,  z^--1)e^-z- (x/2)^2/z z
    )^2   y  .
 
We claim that for <0 we simply have 

    T_1(x) x ↘ 0∼( x/2)^∫_0^∞ y^ -1e^y( ∫_y^∞( z^-  -|| z^--1)e^-z z
    )^2   y  =:  2^||Γ(||)x^   .
 
For this  choose ∈(0,1) and split  the integral in y
in (<ref>) as ∫_x/2^∞…=  ∫_x/2^… +   ∫_^1/… +  ∫_x/2^1/…:= I_1+I_2+I_3. The limit x ↘ 0 is easily taken in I_2 and the dependence on x disappears. Moreover we directly check that lim_↘ 0lim_x↘ 0I_2∈ (0, ∞) is the integral in the right-hand side of (<ref>). 
We are left with showing that lim_↘ 0sup_x∈(0, 2) I_j=0 for j=1 and 3.
For I_1 recall that in (<ref>)  we can replace ∫_y^∞… z with  ∫_0^y … z,  so that 
for  sufficiently small (so x is small too and we can use the asymptotic approximation of _, ∼||) we have 

    sup_x∈(0, 2)| I_1| ≤ ∫_0^ y^ -1( ∫_0^y ( z^-  +2 ||  z^--1)  z
    )^2    y ↘ 0 ⟶ 0 .

Moreover (again,  small)

    sup_x∈(0, 2)| I_3| ≤  2
    ∫_1/^∞ y^ -1e^y( ∫_y^∞( z^-  + 2 || z^--1)e^-z z
    )^2   y ≤
    ∫_1/^∞ e^3y/2( ∫_y^∞ e^-z z  
    )^2   y ↘ 0 ⟶ 0  ,

and (<ref>) is proven.



For =0 we have 

    T_1(x)   = ∫_x/2^∞ y^-1e^y+ (x/2)^2/y( ∫_0^y ( 1  -_, 0 z^-1)e^-z- (x/2)^2/z z
    )^2   y  .
 
We anticipate (for future use) that the result we are going to obtain would be the same if 1  -_, 0 z^-1 is
replaced by _, 0 z^-1
Again,  ∫_0^y… can be replaced by ∫_y^∞…  and it suffices the splitting
∫_x/2^∞… =  ∫_x/2^… +  ∫_^∞… =:I_1+I_2.
In fact (recall that _,0=o(1) as x ↘ 0, in particular _,0 becomes smaller than one) 

    sup_x∈(0, 2)| I_2 | ≤  2∫_ ^∞
     y^-1 e^y (∫_y^∞ (1+ z^-1)e^-z z)^2
       y  ,

and the right-hand side is just a finite expression that depends on .
On the other hand I_1 diverges as x ↘ 0. In fact observe that

    e^-2 Ĩ ≤  I_1  ≤  
    e^2 Ĩ ,  with Ĩ  := ∫_x/2^ y^-1( ∫_0^y ( 1  -_, 0 z^-1)e^- (x/2)^2/z z
    )^2   y  ,

so we can focus on Ĩ. By using ∫_0^L (1/x)exp(-1/x)  x ∼log L for L →∞,
and the fact that /x^2≤ y/x^2≤ 1/(2x) we see that 

    ∫_0^y  z^-1 e^- (x/2)^2/z z ∼log(y/ x^2)  ,

uniformly in the range of y we are using, and as x ↘ 0. Therefore

    Ĩ∼( 1/2 log(1/x))^2 ∫_x/2^ y^-1( 2 log(1/x) - log (1/y))^2  y ∼ 7/12log (1/x) .

This concludes the =0 case: T_1(x) ∼7/12log (1/x).



The case >0 is quicker to treat for >1/2 because of Remark <ref>. In reality also for ∈ (0,1/2] the term containing _, does not contribute: we will check  this fact after  estimating what is giving the main contribution:

    ( x/2)^∫_x/2^∞ y^ -1e^y+ (x/2)^2/y( ∫_y^∞  z^- e^-z- (x/2)^2/z z
    )^2   y  ,
 
and the final result is that for ∈ (0,2)

    T_1(x)  ∼ ( x/2)^∫_0^∞ y^ -1e^y( ∫_y^∞  z^- e^-z z
    )^2   y  =:  q_1() x^ .

This is proven like before by restricting the integral to y ∈ (,1/) and estimating the rest before 
letting ↘ 0. The function q_1(·) can be expressed with a  Meijer G-function, but this does not make it 
much more explicit. 
Let us quickly verify that the term we neglected for ∈ (0, 1/2] is of lower order: 
by focusing on y∈ (x/2, ) (otherwise the fact is obvious) we see that an upper bound on this contribution
is O(x^5) ∫_x/2^ y^- -1 y = O(x^4).




The =2 case generates a logarithmic correction: in fact from (<ref>) we see that if we restrict the integral over y≥, the contribution is bounded by x^2 times a constant that depends only on . 
The integral with y∈ (x/2, ) instead is controlled above and below, up to a factor that can be chosen 
arbitrarily close to one uniformly in x↘ 0 by choosing  small (like in (<ref>)), by

    ( x/2)^2∫_x/2^ y ( ∫_y^∞ z^-2e^-z z)^2  y
    x ↘ 0∼( x/2)^2log (1/x) .




For >2 we go back to (<ref>)

    ( x/2)^∫_x/2^∞ y^ -1e^y+ (x/2)^2/y( ∫_y^∞  z^- e^-z- (x/2)^2/z z
    )^2   y  = ( x/2)^∫_1^∞… + ( x/2)^∫_x/2^1 …≤
     
    C x^ + 
     3( x/2)^∫_x/2^∞ y^ -1( ∫_y^∞  z^- e^-z z
    )^2   y ≤  C' x^(1+ ∫_x/2^∞ y^- +1 y)  =  O(x^2)  ,

where C and C' are constants independent of x.



We collect what we have obtained:

    T_1(x) ∼ 
    2^||Γ(||)
     x^    if ∈ (-∞,0)  ,
    7/12log (1/x)     if =0  , 
    
    q_1() x^    if ∈ (0,2)   ,
    1/4 x^2log (1/x)     if =2  , 
    
    O(x^2)    if >2  .




We now turn to T_2(x) and the basic expression is after a change of variables  (still, =√(2))

    T_2(x)  =
    ( x/2)^4-∫_x/2^∞ u^-1- e^u+(x/2)^2/u(
     ∫_u^∞(  v^-2+ 
     - ( 2 _, /x^2) v^-1+) e^-v-(x/2)^2/v v
      )^2  u  .
 

Let us start with <0 and recall that by Remark <ref> it suffices to consider 

    ( x/2)^4-( 2 _, /x^2) ^2
    ∫_x/2^∞ u^-1- e^u+(x/2)^2/u(
     ∫_u^∞   v^-1+ e^-v-(x/2)^2/v v
      )^2  u  .

The pre-factor behaves asymptotically as ^2 (x/2)^- and the integral can be bounded by two times 

    ∫_x/2^∞ u^-1- e^u(
     ∫_u^∞   v^-1+ e^-v v
      )^2  u ∼ 1/||^3(x/2)^  .

So T_2(x)=O(1) for <0.



For =0 the expression to evaluate is

    _,0^2 ∫_x/2^∞ u^-1 e^u+(x/2)^2/u(
     ∫_u^∞   v^-1 e^-v-(x/2)^2/v v
      )^2  u  .

But this term is minimally different from
(<ref>) (see observation right after(<ref>)) and exactly in the same way we arrive at
T_2(x) ∼ T_1(x) ∼7/12log (1/x). 



For ∈ (0,2) we split the integral with respect to u and the contribution when u≥ 1 is bounded 
so the contribution to T_2(x) is O(x^4-). For u<1 we make an upper on the  contribution to T_2(x):

    3( x/2)^4-∫_x/2^1 u^-1-(
     ∫_0^u (  v^-2+ 
     + 2 _, /x^2  v^-1+) e^-v v
      )^2  u
       
    ≤    C x^4-(
      ∫_x/2^1 u^-3+ u+ ( 2 _, /x^2)^2 ∫_0^1 u^-1+ u)  ≤  C' ( x^2  + _, ^2 x^-)  ,
 
and since _, = max(x^2, x^2), except for a logarithmic correction for =1,  
we conclue that T_2(x)=O(max(x^3, x^2)) for ∈ (0,2).



For =2 we again split the integral with respect to u≥ and u<. The integral for 
y≥ is bounded by a constant that depends only on . 
Arguing as in  (<ref>) we see that what it suffices to control

    ∫_x/2^  u^-3(
     ∫_0^u (  1
     - ( 2 _, 2/x^2) v ) e^-v v
      )^2  u ∼ log (1/x) ,

where we have used that _, 2=O(x^2). Therefore T_2(x) ∼ (x/2)^2  log (1/x) for =2.




Finally, for >2 we have

    T_2(x) ∼( x/2)^4-∫_0^∞ u^-1- e^u(
     ∫_u^∞(  v^-2+ 
     -  v^-1+/-1) e^-v v
      )^2  u =  2^-4Γ(-2)/(-1)^2 x^4-.

The proof of this claim 
follows the same line as the proof of (<ref>), that is, splitting of the y integral in three parts and taking 
the limit ↘ 0. 



We have got to:
 
    T_2(x) ∼ 
     O(1)     if <0  , 
    7/12log (1/x)     if =0  , 
    
     O(max(x^3, x^2))      if ∈ (0,2)  , 
    1/4 x^2 log (1/x)     if  x=2  , 
    2^-4Γ(-2)/(-1)^2 x^4-    if >2  .

 Therefore only T_2 contributes to the final result  for >2. Otherwise only  T_1  contributes, except 
 at =0 and 2 where they both contribute and exactly with the same amount: 
 
    T_1(x) +T_2(x)  ∼ 
     2^||Γ(||) x^    if ∈ (-∞,0)  ,
    7/6log (1/x)     if =0  , 
    
    q_1() x^    if ∈ (0,2)   ,
    1/2 x^2log (1/x)     if =2  , 
    2^-4Γ(-2)/(-1)^2 x^4-    if >2  .

The final
 result, i.e. (<ref>), is recovered by dividing by K_(x) and using
K_(x)∼ x^-||Γ (||)/2^1-|| for ≠ 0 and K_0(x) ∼log (1/x).
The constant
C() in (<ref>) is

    C()  :=
    2     if  <0  ,
    
    
    7/6     if  =0  ,
    
    
    q_1() 2^1-||/Γ(||)     if ∈(0,2)  ,
    
    
    1/4     if  =2  ,
    
    
    1/( 8(-1)^3 (-2) )     if  >2  ,

with q_1()  given in  (<ref>).
The proof of Proposition <ref> is therefore complete.


 

§ LYAPUNOV EXPONENT AND SINGULARITIES: THE PROOF OF PROPOSITION <REF>
  
 
 
In view of (<ref>) we just consider ≥ 0.
We treat first the non integer case.



  §.§.§ The case ∈ (0, ∞) ∖

By the connection formula with the other modified Bessel function 
I_ (x), 
we have 10.274 and 10.252

    K_ (x)  = π/2 sin(π)(I_-(x)-I_(x))  ,

with

    I_(x)  := ( x/2)^∑_k=0^∞(  x^2/4)^k/k!  Γ(+k+1)  =:  ( x/2)^Ĩ_(x)  ,

where Ĩ_(x) is a non standard notation, but it singles out the analytic part of 
the I_(·): in fact Ĩ_(·) is an entire function. 
By elementary manipulations we  obtain

    x K_ -1(x)/K_(x)  =  2 
    (x/2)^2( Ĩ_-1(x)/ Ĩ_-(x))-
    (x/2)^2( Ĩ_-+1(x)/ Ĩ_-(x))
    /1-(x/2)^2( Ĩ_(x)/ Ĩ_-(x))  .

Therefore, aiming at expanding this expression for x ↘ 0 up to the first singular term, we obtain

    x K_ -1(x)/K_(x)    = 
    -2 (x/2)^2( Ĩ_-+1(x)/Ĩ_-(x))+
    (x/2)^22Γ(1-)/Γ() 
    +O(x^4)+O(x^2 +2)
    
       =  p_, ⌊⌋(x^2)+ (x/2)^22Γ(1-)/Γ() 
    + O(x^ 2⌊⌋ +2)+O(x^4)  ,

where p_, j(y) is the Taylor expansion up to degree j of 
p_(y):=- (y/2) Ĩ_-+1(√(y))/ Ĩ_-(√(y)). 
It is not difficult to realize that 
the coefficients of this Taylor expansion are  just rational function of .
Let us detail this point that is going to be important also for the passage to  integer:
if we introduce for k∈∪{0} the Pochhammer's symbol 

    (ν)_k  := Γ(ν+k)/Γ(ν)k=1,2, …= (ν+k-1) (ν+k-2) ⋯ (ν+1) ν  ,

we see that p_(y) can be written in terms of Pochhammer's symbols:

    p_(y)  =  - y/2(∑_k=0^∞y^k/k!  (-+1)_k+12^2k)/
    (∑_k=0^∞y^k/k!  (-+1)_k2^2k)
    .
 
From now the explicit determination of p_(y) is elementary, but cumbersome (to the point of requiring symbolic 
computations).
 We give the first four terms

    p_(y)  =    ∑_j=1^4 c_j()y^j + 
     … = 1/2( -1) y
    - 1/8(-2)(-1)^2y^2 +
         1/16  (-3)(-2)(-1)^3 y^3 
    -5  -11/128 ( -4)(-3)(-2)(-1)^4 y^4
    
    + …

This completes the proof for the non integer case.



  §.§.§ The case =0,1, 2, …

We need to treat separately the case =0 because it involves K_-1(x), which
however it is just K_1(x), but it requires an ad hoc (much simpler) analysis. So we start off
with the case =1, 2, …. From now  will be replaced by n and when we write  we mean a 
quantity that is not integer.
 By  10.311 we have that for n=0,1,2,…
 
    ( x/2)^n K_n(x)  = 1/2∑_k=0^n-1 (-1)^k
     (n-k-1)!/k!( x/2)^2k+ (-1)^n/n!( x/2)^2nlog (x) 
     +O(x^2n)  ,

 and the sum as to be interpreted as empty if n=0.
 For n=1, 2, … we write the degree 2n-2 polynomial in the right-hand side as a n-1 degree polynomial 
 with argument x^2:
 
    q_n(x^2)  := 1/2∑_k=0^n-1 (-1)^k
     (n-k-1)!/k!  2^2k( x^2)^k .

 With this notation we have for n=1, 2, …
 
    x K_n-1(x)/K_n(x)    = 
     2 (x/2)^2 q_n-1(x^2)/q_n(x^2) + 2(-1)^n/(n-1)!(x/2)^2n/q_n(x^2)log x + O( x^2n)
     /1+O( x^2n|log x|)
     
        = x^2/2 t_n(x^2)+
     2^2-2n(-1)^n/((n-1)!)^2x^2nlog x + O( x^2n)  ,

 where t_1(y):=0 and, for n=2, 3, …, t_n (y) is the polynomial of degree n-2 
 given by the Taylor expansion of the rational function q_n-1(y)/q_n(y).  
 
 We are therefore left with showing that the coefficients of the polynomial y t_n(y)/2 of degree n-1 coincide with the corresponding coefficients of the Taylor polynomial of
 p_(y), when =n (we can consider n=2,3,… becaue if n=1 the polynomial is identically zero and our claim  is trivially verified).
 In different terms, we have to show that 
 if we set =n=2,3 , …
 in (<ref>) up to the degree n-1 (it is readily seen that
 the n-th Taylor coefficients diverges as → n), then 
 we obtain y t_n(y)/2. To prove this
 it is useful to remark that, in order to obtain the Taylor expansion of p_(y)
 up to order ⌊⌋ it is sufficient to expand the rational function
 
    p̃_(y)  =  - y/2(∑_k=0^⌊⌋-1y^k/k!  (-+1)_k+12^2k)/
    (∑_k=0^⌊⌋-1y^k/k!  (-+1)_k2^2k)
    .
 
In this expression we can set  = n obtaining thus the rational function

    p̃_n(y)  :=  - y/2(∑_k=0^n-2y^k/k!  (-n+1)_k+12^2k)/
    (∑_k=0^n-2y^k/k!  (-n+1)_k2^2k)
    .
 
and Taylor coefficients up to degree n-1 of this function are precisely the limit for → n
of the Taylor coefficients up to degree n-1 of p_(y), cf. (<ref>)-(<ref>).
We are left with showing that the coefficients up to degree n-1 of p̃_n(y)
coincide with the corresponding coefficients of y t_n(y)/2. This is equivalent to showing that 

    -
    (∑_k=0^n-2y^k/k!  (-n+1)_k+12^2k)/
    (∑_k=0^n-2y^k/k!  (-n+1)_k2^2k)  = q_n-1(y)/q_n(y)
    + O ( y^n-1)

and this is implied by the stronger (non asymptotic) condition

    -
    (∑_k=0^n-2y^k/k!  (-n+1)_k+12^2k)/
    (∑_k=0^n-2y^k/k!  (-n+1)_k2^2k)  =
    (
    ∑_k=0^n-2 (-1)^k
    (n-k-2)!/k!   2^2k  y^k)
     /
     (
    ∑_k=0^n-2 (-1)^k
    (n-k-1)!/k!   2^2k  y^k)
     
       ,

For n=2, 3, … this identity is verified directly by using that for n=1,2, … and k∈ 0,1, …

    (-n)_k =  (-1)^k n!/(n-k)!  ,

and this completes the proof in the case of =1,2, ….

We are left with the case =0:

    x K_-1(x)/K_0(x)  = x K_1(x)/K_0(x)(<ref>)∼1/log (1/x) ,

but this is easily improved by going back to 10.311 and using

    K_0(x)  =  -log (x/2) - γ + O(x^2)  ,

 where γ is the Euler-Mascheroni constant. This, with (<ref>) for n=1, implies that
x K_-1(x)/K_0(x) is equal to 1/(log(1/x) +(log 2- γ))+ O(x^2).
This completes the proof of Proposition <ref>.


 

 

§ LYAPUNOV EXPONENT AND SINGULARITIES: THE PROOF OF THEOREM <REF>

 
 
 The proof of Theorem <ref> is somewhat involved, since the ratio of Bessel functions in the integrand becomes quite singular at one end of the domain of integration. 
 For every fixed x>0 the numerator and denominator are entire functions of , so the ratio is analytic apart from the zeros of the denominator; 
 these are all on the imaginary axis and they are bounded 
 away from the real axis as long as x is bounded away from zero (<cit.>, Table <ref>).
 But integrating over x down to zero adds contributions that are
 less and less regular as the gap between the origin and the poles of the integrand 
 shrinks when x becomes small. 
  
 
 
One way to obtain a proof is to exploit  
once again  the
  connection formula (<ref>)-(<ref>) which gives an expansion 
 of both numerator and denominator:  but controlling the ratio is of course not straightforward. 
 And in fact  McCoy and Wu approach the problem this way, but keeping only the 
leading terms of the series in the connection formula. 
The validity of this procedure is not obvious, since a priori the resulting correction could be less regular than the leading terms, 
but it is nonetheless helpful to begin by examining this simplified problem
that has the nice feature of leading to a solution in terms of special functions, since we shall see that it correctly illustrates the main features of the proof. 
With this aim in mind we examine
a simplified McCoy-Wu formula – this corresponds to studying a function
 which is defined, like , as the integral over x∈ (0, η) of a suitable function
f̃_x() (see (<ref>)) – and



  * we perform the integration explicitly and discuss the regularity  of ↦(); 

  * we then argue how understanding the location of (some of) the poles
of f̃_x in the complex plane, and the corresponding residues, gives another way to understand the regularity.
  
 
 All of this is done in Section <ref>. Then 
in Section <ref> we give   the proof of Theorem <ref>, based on a treatment of the poles of f_x.
 


 §.§ Heuristic arguments and idea of the proof




  §.§.§ The simplified McCoy and Wu problem: exact solution

Much like McCoy and Wu did in  <cit.>, we can consider the  
leading contribution to the integrand in (<ref>) for
x↘ 0 and ∈ tending also to zero. 
This step is an uncontrolled approximation that is obtained by keeping 
the first terms in (<ref>)-(<ref>)
and by using Γ() ∼ 1/ 

    K_ - 1(x)  ≃ 1/2( x/2)^ -1  ,

and 

    K_(x)  ≃ 1/2( ( x/2)^- - ( x/2)^)  ,

 so
 
    x K_-1(x) /K_(x) ≃2 /( 2/x)^2 -1  = 2 /exp(2  L(x))-1  =:  f̃_x(α)   ,

 with L(x):= log(2/x). 
 We then have
 
    ()  := ∫_0^η2 /exp(2  L(x))-1 x  .

 for η∈ (0,2), cf. <cit.>.
 From now onward the analysis is rigorous.
 
 
 
 
 
() can then be made explicit  up to an additive contribution that is analytic near the real axis: in fact
if we set

    ()  := ∫_0^2  (2 /exp(2  L(x))-1 -1/L(x)) x  ,

we have

    ()  = ()
     -
      ∫_η^2  (2 /exp(2  L(x))-1 -1/L(x)) x +
    ∫_0^η1/L(x) x   .

But ∫_0^η (1/L(x))  x
= Γ(0, log(2/η)) is just a constant
and the second addend in the right-hand side is (real) analytic in . In fact
 the integrand in  is meromorphic with poles on the imaginary axis, precisely for  equal to any integer multiple of ±π /L(x). Therefore the second addend in the right-hand side is  analytic for 
in ∖{iy: | y|≥π /L(η)}. We can therefore focus on  
() and we start by observing that for ∈ 

    (-)  =  ()+ 4 ,

 as a result of a straightforward manipulation using 
 1/(1-e^-2 L)= 1/(e^2 L -1)+1.
 
 
 
 We note that the same argument can be applied directly to _1 () obtaining
 
    (-)  = ()+ 2 η .

 But in fact we have also 
 
    (-)  =  ()+ 2 η ,

 which follows by applying  (<ref>) and  K_(x)=K_-(x). 
 In other words all these functions – ,  and   – are even, up to a linear term.
 
 
 Thanks to (<ref>) we  can focus on the case  >0
 and compute (change the variable and move the contour of integration in the complex plane)

    ()  =  2 ∫_0^∞(1/e^v -1- 1/v) e^-v/(2) v  =  
     -2 ∫_0^∞( 1/v - 1/1-e^-v) e^-v/(2) v -4 
    
    =  -4  -2 log (2 ) -2 ψ(1/(2)) α↘ 0∼  -2  - ∑_j=1^∞B_2j/j (2 )^2j   ,

where ψ(z)= Γ'(z)/Γ(z) and the notion of ∼ is extended here
and it has to be interpreted in the sense of asymptotic series (i.e., that the difference of left-hand side and of the series in  the right-hand side
truncated to j=n is o(^2n))
: in the third step we have applied 5.913  

and the asymptotic relation  5.112 
ψ(z)z → +∞∼log z -1/(2z) + ∑_j=1^∞ (B_2j/(2j))z^2j: the rational numbers B_2n are the Bernoulli numbers 24.2. By 24.98	

    B_2nn →∞∼ 2 (-1)^n+1(2n)!/(2π)^2n  ,

so the series has radius of convergence zero. Note that
 (<ref>) readily implies that
()∼ -2  - ∑_j=1^∞ (B_2j/j) (2 )^2j 
holds also for ↗ 0 and not only for ↘ 0. 
 Let us reorder what we have done (and more) into a statement:





(·) is defined and analytic in the complex plane without the imaginary axis
and it can be extended by continuity on the whole real axis  by setting (0)=0. Then
 restricted to  is C^∞ in the origin. On the other hand,  cannot be 
continued as an analytic function at any point on the imaginary axis.




Proof.
Let us first show that
  is C^∞ at the origin. For this it is more practical  (and equivalent, since -  is real analytic) to go back to 
  and observe that with y=2∈∖{0}
 and the change of variable v=log (2/x) we obtain
 
    f(y) := (y/2) =    2∫_c^∞ q(yv) exp(-v)/v v  ,
 
 with  c=log (2/η)∈ (0, ∞)  and 
 q(u) = u/exp(u)-1.
 It is straightforward to verify that q(u)∈ (0,1] for u≥ 0 and that q(u) ∈ [1, u+1] for u≤ 0
 and dominated convergence implies that f is C^0 (on the whole , but of course our attention 
 is at the origin, outside we already know that f is  real analytic).
 To show differentiability it suffices to observe that formally 
  
    f^(n)(y) = ∫_c^∞ q^(n)(yv) v^n-1exp(-v) v  ,
 
 where  f^(n)(y)= ( / y)^n f(y). But  it is straightforward to verify that q^(1)(u) is monotonically decreasing from 
 1 to zero and q^(n)(u), n=2 or larger, vanishes at ±∞. Hence sup _u ∈| q^(n)(u)|
 < ∞ and this implies that 
 (<ref>) is not just formal: f ∈ C^∞ and its derivatives are given by (<ref>).
 
 
 
 
 On the other hand,  is  not analytic in zero since the radius of convergence of the series is zero: 
we record that if the Taylor series for () is
∑_n c_n  ^n then  

    c_2n n →∞∼  4 (-1)^n+1(2n-1)!/(π)^2n  ,

while c_2n+1=0 for n=1,2, …: this follows directly from (<ref>) and (<ref>).

To conclude, we argue that  cannot be extended as an analytic function through the imaginary axis. 
Call ^+: {z ∈:  ( z) >0}⟶ the function that coincides  with  in its domain of definition (the right half-plane): keep in mind that  is defined also in the half plave with negative real part and it is analytic there.  
Suppose now that there exist a∈ and ϵ >0 such that we can define ( ) for ∈{it:  | t-a | <ϵ}, in such a way that  is analytic in a neighborhood of ia. We can and do assume that a>0 as well as a->0 by symmetry and because we already know that 
is not analytic at the origin. 
This extends  analytically ^+ too. 
But (<ref>) yields ^+() = -4  - log 4-2 log -2 ψ(1/(2)) in the right-half plane and this expression, containing 
the function  ψ (that is meromorphic on the whole  and has poles on the negative real axis 5.2) and a logarithm which is defined and analytic 
on the whole complex plane except for a cut (that can be for example chosen to be {z∈:  z = z ≤ 0}).
But ^+ must coincide with  in the whole region where ^+ is extended: this region includes the
negative semi-axis where  ^+ has poles and  is analytic. This is not possible,
 hence  cannot be extended 
analytically at any point on the imaginary axis.








  §.§.§ A different viewpoint on the McCoy and Wu simplified problem

We restart from (<ref>), but we take a different approach: we avoid exact integration. 
As already noted,
the integrand in that expression for (), for fixed x>0,
is a meromorphic function of . The poles 
are on the imaginary axis and 0 is not a pole because the singularity is removable:
the poles are n π  i / L(x) for n ∈∖{0}, where recall that we write L(x) := log (2/x) for brevity.
The  integrand in the expression for _1 () is therefore analytic in a ball of radius   π / L(x) around zero
and the most singular part in the residue expression for the integrand comes from the two closest poles, which are 
±π i/ L(x). If we compute the residues of these two poles for the integrand  we find ± 2π i /(L(x))^2 and therefore
the contribution to the integrand of these two poles is 

    P_(x)  := 2π i/L(x)^2 (2 - 2π i/L(x))- 2π i/L(x)^2 (2 + 2π i/L(x))    = 
    - 2/L(x)1/( L(x)/π)^2+1
        =  -2 ∑_n=0^∞ (-1)^n ^2n L(x)^2n-1π^-2n  ,

where the last equality holds only for | L(x)/ π| < 1, but it is in any case useful to identify 
all the derivatives of P_ (x) at =0. It seems reasonable to believe that the singularity at the origin
of  is induced by the poles of the integrand and that the two poles that are closest to the origin
give the leading part of the singularity. If this is the case ↦∫_0^η p_(x)  x should capture
the leading behavior of the singularity of . This is confirmed or at least highly suggested 
by observing that

    ∫_0^η L(x)^2n-1 x   =  2 ∫_log (2/ η)^∞ y^2n-1 e^-y y   =  2 Γ(2n, log (2/η))
    n →∞∼  2 Γ(2n)  =  2 (2n-1)!  ,

so if we proceed at a completely 
 formal level, by  integrating term by term the series in the second line of   (<ref>) and using (<ref>),
 we directly recover 
 (<ref>)!
 



  §.§.§ Approaching the true problem


Going back to the true problem, that is , for which the integrand is  the left-had side of (<ref>), we have to identify the zeros 
of K_(x): this problem has been studied and, as possibly expected, the heuristics coming from studying the poles 
of  the right-hand side of (<ref>) is qualitatively correct, so the poles we have to study, or the zeros of  
K_(x), are on the imaginary axis. Moreover they accumulate on the origin when x↘ 0. In fact, to leading order 
(for n fixed and x ↘ 0) they are still in ± n π i  /log (2/x). But, more precisely,
they are in ± n π i  /(log (2/x)- γ+ o(1)) as x ↘ 0 (γ the Euler constant, see below). 





The subleading correction -γ to log (2/x) for the location of the zeros can be inserted in the 
heuristic argument that we presented just by being keeping one more term in the expansion of the Γ function in the first steps 
(<ref>)-(<ref>).
The change amounts simply 
 to work with L(x)= log (2e^-γ/ x), but this small offset in L(x) leads to the multiplicative 
e^-γ constant in the final asymptotic results, see Theorem <ref>. 



Apart for the quantitative issue of Remark <ref>, the proof
of Theorem <ref> requires taking care of two main issues:





  * Even admitting that the two closest poles give the main contribution, we have to set-up a rigorous procedure
corresponding to the formal argument that we developed using (<ref>) and (<ref>).
This procedure requires  controlling not only location of the two poles, but also the residues that
this time are given by ratio of series coming from (<ref>).

  * One needs to control the effect of the poles that are farther from the origin: note that these terms will in any case give contributions
that  generate divergent series, but this time the distance of the poles from the origin is at least (about) twice the distance of the two n=1
poles. Hence they will contain an exponential factor that is at least twice smaller. 








 §.§ The proof
 


The  proof of Theorem <ref>, that starts here, follows a main line separated by a number of lemmas and corollaries.   
As a preliminary, we set out precise versions of some of the statements made above about the regularity of 

    f_x: α↦ xK_1-α(x)/K_α(x),

for x > 0.
K_α(x) is entire as a function of α for all x ≠ 0 <cit.>,
so f_x is a ratio of two entire functions, and thus it is analytic except at the zeros of α↦ K_α (x).  For x>0, these zeros are all pure imaginary and located in the region |α| > x <cit.>.
Note that

this characterization
is already sufficient to show that ∫_X^η f_x (α) d x is analytic on ∖± i [X,∞) for any 0 <X < η < ∞, so it suffices to prove the statements in Theorem <ref> with η replaced by some sufficiently small X in the integral defining . 
Furthermore, f_x is infinitely differentiable on the real axis for any x > 0; in the next few lemmata we will bound its derivatives in such a way as to show that this is also true of the integral .

To begin with, note that combining (<ref>) and (<ref>) gives

    K_α(x)
    	=
    	π/2 sinπα∑_k=0^∞( x^2 / 4 )^k/k![ 
    		( x/2 )^-α/Γ(-α+k+1) 
    		- ( x/2 )^α/Γ(α+k+1)],

for α∈∖ and x > 0.



	There exist A,C_1,X > 0 such that for all x ∈ [0,X], α∈ such that |α|≤12 and |α| log 1/x ≤ A,
	
    |x^α/Γ(1+α) - x^-α/Γ(1-α)|
    		≥ C_1 |α|x^-|α|.

	



	Letting g_x(α) := x^α/Γ(1+α),
	
    x^α/Γ(1+α) - x^-α/Γ(1-α)
    		=
    		g_x(α)-g_x(-α)
    		=
    		α∫_-1^1 g'_x(yα)  y   ,

	and so
	
    |x^α/Γ(1+α) - x^-α/Γ(1-α)|
    		≥ 
    		|α|
    		|
    		∫_-1^1  g'_x(yα)  y
    		|  
    		.

	Note that
	
    g'_x(α)
    		=
    		[ log x - ψ(1+α) ]
    		x^α/Γ(1+α)
    		,

	where ψ(z) := Γ'(z)/Γ(z) is the Psi function 5.2.
	ψ(1+α)/Γ(1+α) is analytic for |α| < 1, so choosing X small enough we can obtain | ψ(1+α)/Γ(1+α) | ≤ (1-c_1) for any fixed c_1 ∈ (0,1), whence
	
    | g'_x(α) - x^αlog x /Γ(1+α)|
    		≤|g'_x(α) - x^αlog x /Γ(1+α)|
    		=
    		| ψ(1+α) x^α/Γ(1+α)|
    		≤
    		(1-c_1) 
    		x^αlog 1/x
    		,

	and noting that 1/Γ(1+α) is entire and takes positive real values for α∈ [-1/2,2/2], by choosing A small enough we obtain
	
    -  x^αlog x/Γ(1+α)
    		=
    		x^αlog 1/x [ cos( αlog x ) 1/Γ(1+α) 
    			- sin( αlog x ) 1/Γ(1+α)]
    			
    ≥ c_2 x^αlog 1/x ≥ 0  ,

	for some c_2>0.
	Combining the above observations, we see that
	
    |
    		∫_-1^1  g'_x(yα)  y
    		|
    		≥ 
    		c_1 c_2
    		log 1/x
    		∫_-1^1 x^y α y
    		=
    		c_1 c_2
    		x^-α - x^α/α.

	The last expression is an even function of α, so without loss of generality we can consider α = ρ∈ [0,1/2].  Noting that 
	ρ↦ 1 - x^2ρ = 1- exp(- 2 ρlog x ) is concave for all x > 0 and checking the boundary cases, it is easy to see that 1 - x^2ρ≥ρ for all x ∈ (0,1/2), ρ∈ [0,1/2].
	Using this together with (<ref>)  and (<ref>), we obtain
	
    |x^α/Γ(1+α) - x^-α/Γ(1-α)|
    		≥ 
    		c_1 c_2
    		|α| x^-|α|.

This completes the proof of Lemma <ref>




	There exist some X ∈ (0,1), A > 0, and C_2 > 0 such that 
	
    |f_α (x)| ≤  C_2  ,

for all x ∈ (0,X], |α| ≤ 1/2, and |α| ≤ A/log(2/x).
	


	Rearranging (<ref>), we have
	
    - 2 sinπα/πK_1-α(x)
    			   =
    			(x/2)^α-1/Γ(α)
    			+
    			∑_k=1^∞( x/2)^α+2k-1/k!Γ(α+k)
    			-
    			∑_k=1^∞k/k-α( x/2)^-α+2k-1/k!Γ(k-α)
        =
    			(x/2)^α-1/Γ(α)
    			+
    			∑_k=1^∞(x/2)^2k-1/k![(x/2)^α/Γ(k+α) - k/k-α(x/2)^-α/Γ(k-α)],

	Noting that the final sum vanishes term by term when α=0 and that
	
    ∂/∂α[ (x/2)^α/Γ(k+α)]
    		=
    		(x/2)^α/Γ(k+α)[ log (x/2) - ψ(k+α) ]  ,

	and
	
    ∂/∂α[ (x/2)^-α/Γ(k-α)k/k-α]
    		=
    		 ∂/∂α[ k (x/2)^-α/Γ(1+k-α)]
    		=-
    		k (x/2)^-α/Γ(1+k-α)[ log(x/2) + ψ(1+k-α) ],

	and that |ψ(k+α)|∨|ψ(1+k-α)| ≤log (k+1) for k ≥ 1 and |α| ≤ 1/2, we then have also
	
    | (x/2)^α/Γ(k+α) - k/k-α(x/2)^-α/Γ(k-α)|  ≤
    |α| [ |(x/2)^-α/Γ(k+α)| ∨|k (x/2)^α/Γ(k-α)| ]
    		[ log(x/2) + log (k+1) ].

	Noting that Γ(1+k ±α) = Γ(1 ±α) (1 ±α)_k (see (<ref>)),
	
    |(1 ±α)_k|
    		=
    		|1 ±α| ⋯|k ±α|
    		≥( 1/2) ⋯( 2k-1/2)
    		= 2^-k  (2k-1)!!,

	and that |Γ(1 ±α)| is bounded for |α|≤ 1/2, we have 
	
    | ∑_k=1^∞(x/2)^2k-1/k![(x/2)^α/Γ(k+α) - k/k-α(x/2)^-α/Γ(k-α)] |
    			
         ≤ |α| ∑_k=1^∞(x/2)^2k-1/k![ log(k+1) + log (2/x) ] [ | (x/2)^-α/Γ(k+α)|∨|  k (x/2)^α/Γ(k-α)| ]
    			
         ≤ |α| (x/2)^1-|α|[ ∑_k=0^∞2^k+1log (k+2)/k!(2k+1)!!(x/2)^2k 
    			+ log(2/x) ∑_k=0^∞2^k+1(x/2)^2k/k!(2k+1)!!]
    			
         ≤ |α| [ 1 + log (2/x) ] (x/2)^1-|α|,

	where the last bound follows from the observation that each of the sums in the preceeding expression is α-independent and defines an entire function of x, and is therefore bounded on any compact interval.  Combining this with (<ref>) and noting that
	
    | (x/2)^α-1/Γ(α)|
    		≤ |α| (x/2)^α -1,

	we have
	
    | 2 sinπα/πK_1-α(x) |
    		≤ |α| (x/2)^α -1.


	As for the denominator, using (<ref>)
	
    |2 sinπα/πK_α(x)|
    			≥   | ( x/2 )^α/Γ(1+α) - ( x/2 )^-α/Γ(1-α)|
    			
        -
    			∑_k=1^∞(x/2)^2k/k!| ( x/2 )^α/Γ(k+1+α)
    			-
    			( x/2 )^-α/Γ(k+1-α)|  .

	The sum can be estimated in the same way as the one in (<ref>): we have
	
    | ( x/2 )^α/Γ(k+1+α)
    			-
    			( x/2 )^-α/Γ(k+1-α)|
    			≤2^k [ log(k+1) + log (2/x) ]/(2k-1)!! |α| ( x/2)^-|α|  ,

	and so
	
    ∑_k=1^∞(x/2)^2k/k!| ( x/2 )^α/Γ(k+1+α)
    			-
    			( x/2 )^-α/Γ(k+1-α)|
    			≤ |α| x^2 - |α|log (1/x)  ,

	and we see that this is dominated by the first term, which was bounded from below in Lemma <ref> above.  We then have
	
    |2 sinπα/πK_α(x)|
    		≥× |α| ( x/2 )^-| α|,

	and combining this with (<ref>) we obtain the desired bound
	and the proof of Lemma <ref> is complete.


Letting C_R(w) denote the oriented circle of radius R about w, the Cauchy formula implies that 

    |f^(n)(w)|
    	=
    	n!/2π| ∮_C_R(w)f(z)/( z-w )^n+1 d z|
    	≤n!/R^nmax_|z-w|=R|f(z)|  ,

for any f which is analytic on an open set containing C_R(w) and its interior.  From Lemma <ref> we thus have

 
	For the same A,X,C_2 as in Lemma <ref>, 
	
    |∂^n/∂ a ^n[ x K_1-a(x)/K_a(x)]|
    		≤ C_2 n! ( log (2/x)/A)^n

	for all n ∈, a ∈ [-1/4,1/4].


Since the bounds in Corollary <ref> are uniformly (in a) integrable (in x), if we take η≤ X this allows us to take derivatives inside the integral defining  which is therefore infinitely differentiable on the real interval (-1/4,1/4). This result
is going to be crucial for us at 0: the fact that  is C^∞ outside of zero is also a byproduct  of the fact that
we are going to establish (via the next lemma) that  is real analytic in (-1,1)∖{0}.  





	For any A ∈ (0,1/2) and I ∈ (0,∞), there exist X ∈ (0,2) and C ∈ (0,∞) such that
	
    |f_α (x) |  ≤   C  ,

whenever x ∈ (0,X], α∈ (A,1-A), and |α| ≤ I.
	


	Using 5.67 and noting that Γ(x) > 1/2 for all x > 0,
	
    |1/Γ(z)| 
    		≤( coshπ z )^1/2/Γ( z)≤ 2 ( coshπ z )^1/2  ,

	whenever z > 0.  Also using 5.66 and noting that the Gamma function is concave for positive real arguments,
	
    |1/Γ(2-α)|
    		≥1/Γ(2 - α)≥1/Γ(1)∨1/Γ(2)
    		=
    		1.

	Applying these two bounds to (<ref>), we obtain
	
    |πsinπα/π K_α(x)|
    			   =
    			| 
    			(x/2)^-α/Γ(1-α) 
    			+ ( x/2)^2-α∑_k=0^∞(x/2)^2k/(k+1)! Γ(k+2-α)
       movemovemovemovemove
    			- ( x/2)^α∑_k=0^∞(x/2)^2k/k! Γ(k+1+α)|
    			
       ≥
    			|α| (x/2)^-α/|Γ(2-α)|
    			- 
    			2 [ ( x/2)^2-α +  ( x/2)^α]
    			( coshπ I )^1/2∑_k=0^∞(x/2)^2k/k!
       ≥
    			A ( X/2)^-A
    			-
    			4 ( X/2)^A exp( x^2/4) ( coshπ I )^1/2  ,

	for all relevant α and x; choosing X sufficiently small, the last bound can be made positive.
	

	Similarly, noting also that Γ(x) is negative (resp. positive) and decreasing for x ∈ (-1/2,0) (resp. (0,1/2)), 
	
    |πsinπα/π K_1-α(x)|
    			=
    			| 
    			(x/2)^α-1/Γ(α) 
    			-
    			(x/2)^1-α/Γ(-α) 
    			+ ( x/2)^α+1∑_k=0^∞(x/2)^2k/(k+1)! Γ(k+1+α)
    movemovemovemovemovemove
    			- ( x/2)^3-α∑_k=0^∞(x/2)^2k/(k+1)! Γ(k+1-α)|
    			
    ≤
    			2 ( coshπ I )^1/2{(x/2)^A-1/Γ(A) 
    				-
    				(x/2)^A/Γ(-A) 
    				+ 
    				[ ( x/2)^A+1 +  ( x/2)^A+2]
    				exp(x^2/4)
    			} ,

	and combining this with (<ref>) we obtain a suitable bound on |f_α(x)|.
	This completes the proof of Lemma <ref>.


By Lemma <ref> we have 
that for any α∈ with α∈ (0,1), we can choose A,I to obtain such a bound on a neighborhood of α; using the Cauchy formula this impies that f'_x is uniformly bounded on some smaller neighborhood of α, which allows us to exchange differentiation and integration to see that  is holomorphic on that neighborhood.
We can then conclude that  is analytic on {α∈| α∈ (0,1) }, and by the symmetry noted in (<ref>) it is also analytic on {α∈| α∈ (-1,0) }.



All that remains is to show that the derivatives of  at the origin grow as stated; since this will imply that the associated Taylor series is divergent, this will also prove that  is not analytic there.
We begin by providing a more precise characterization of the poles of f_x for small x.


	
  There exist X,C>0 and a sequence of functions ν_n : (0,∞) ↦ (0,∞), satisfying

    | ν_n(x)  - n π/log(2/x) - γ|
    	≤C n^3 /(log x )^4

and ν_1(x) < ν_2(x) < …,
such that for all x ∈ (0,X], K_α (x) = 0 iff α = ± iν_n(x).


	For x> 0 and ν∈, (<ref>)  can be rewritten using some properties of the Gamma function (5.25 and 5.43) as

    K_iν(x)
    	=
    	-( π/νsinh( πν))^1/2∑_k=0^∞( x^2/4 )^k/k!sin( θ_k (ν) )/√((1^2+ν^2)…(k^2+ν^2))
    	,

where

    θ_k (ν) :=  νlog(x/2) -  Γ(1+k+iν)  ,

(cf. <cit.>, where this expression is used to study the x-zeroes and their dependence on ν).


Then the solutions of K_iν(x) = 0 are the nonzero solutions of

    sinθ_0(ν)
    	=
    	S_x(ν)
    	:=
    	-
    	∑_k=1^∞( x^2/4 )^k/k!sin( θ_k(ν) )/√((1^2+ν^2)…(k^2+ν^2)).



Using the definition of the ψ function and the expansion 5.76,

    θ'_0(ν)
      =
      log(x/2)
      -
      ψ (1+iν)
      =
      log(x/2)
      + γ - ∑_m=1^∞ν^2/m(m^2+ν^2)≤log(x/2)
      +γ,

so for 0 < x < 2e^-γ = 1.1229… θ_0 is strictly decreasing.
From its definition in (<ref>), it is apparent that S_x can be bounded 

    |
    	S_x(ν)
    	|
    	≤∑_k=1^∞( x^2/4 )^k/(k!)^2≤∑_k=1^∞x^2k/(2k)!≤cosh x - 1,

uniformly in ν, and so
we see that for each  0< x < cosh^-1(2) = 1.317…
all solutions ν of 
(<ref>) satisfy

    -θ_0(ν)
    	∈
    	[n π- sin^-1( cosh x - 1  ),n π+ sin^-1( cosh x - 1  )].

for some integer n, and there is at least one solution for each n.
More precisely, for such solutions the derivative of the left hand side of (<ref>) satisfies

    (-1)^n+1∂/∂νsin( θ_0(ν) )
    	=
    	(-1)^n+1θ'_0(ν) cos( θ_0(ν) )
    	≥
    	|θ'_0(ν)|√(1-(cosh x - 1)^2),

and in light of (<ref>), for X small enough this can be bounded from below by any positive number uniformly in ν. As for the right-hand side of (<ref>), first note that from 5.52 we have

    ψ(1+k+α) 
    	= 
    	ψ(1+α) 
    	+
    	∑_m=1^k 1/m+α  ,

for all k ∈, α∈,
and so

    |θ'_k(ν) - θ'_0(ν)|
    		   =
    		|ψ(k+1+α)- ψ(1+α)|
    		≤|ψ(k+1+α)- ψ(1+α)|
    		
       ≤∑_m=1^k 1/m≤∫_1^k+1dm/m
    		= log (k+1).

Using this, we have

    | ∂/∂νsin( θ_k(ν) )/√((1^2+ν^2)…(k^2+ν^2))|
    		=
    		|θ'_k(ν) cosθ_k(ν) - ∑_m=1^k ν/m^2 + ν^2sinθ_k(ν)/√((1^2+ν^2)…(k^2+ν^2))|
    		
    ≤|θ'_0(ν)| + 2 log(k+1)/k!≤|θ'_0(ν)| + 2/(k-1)!,

which implies  

    |S'_x(ν)|
    	≤∑_k=1^∞(x^2/4)^2k/k!|θ'_0(ν)| + 2/(k-1)!≤[  |θ'_0(ν)| + 2]
    	∑_k=1^∞x^2k/(2k-2)!≤[  |θ'_0(ν)| + 2]
    	x^2 cosh x,

which, for x small enough, is smaller than the right-hand side of (<ref>) for all ν.
This suffices to show that there is only one solution of (<ref>) in each of the intervals in (<ref>).  The solution for n=0 must be the trivial solution ν=0 which does not correspond to a solution of K_i ν(x) = 0.
This uniqueness also implies that the solutions for negative and positive n are related by the symmetry K_iν(x) = K_-iν(x), so we see that it is possible to relate the zeros to a family of functions as desired.

To see that the functions ν_n satisfy the bound (<ref>), 
we first note that (<ref>) and (<ref>) together imply that

    ν_n (x) 
      ≤n π + cosh x - 1/inf_ν≥ 0 |θ'_0(ν)| ≤n π + cosh x - 1/log (2/x) - γ≤n/log 1/x  ,

for all x ∈ (0,X);
and also that using the same expansion as in (<ref>) we have

    |
    	γ + ψ(1+iν)
      |
      ≤∑_m=1^∞ν^2/m(m^2+ν^2)≤ν^2 ζ(3)
    	|
    		γν + Γ (1+i ν) 
    	|
    	≤ |ν|^3 ,

for all ν∈, where ζ(s) := ∑_m=1^∞ m^-s is the Riemann zeta function <cit.>.
Then recalling the definition of θ_0, this imples

    |
      [ log (x/2) - γ] ν_n(x)  - θ_0 (ν_n(x)) 
      |
      =
    	|
    		γν_n(x) + Γ (1+i ν_n(x)) 
    	|
    	≤n^3/(log 1/x)^3  ,

and restating (<ref>) as |θ_0( ν_n (x)  ) - nπ| ≤ x^2 this gives the desired bound and the proof of Lemma <ref> is complete.


We denote the residue of f_x at ± i ν_n(x) by ± R_n(x).
Letting 

    K̃_α(x)
    	:=
    	π/2 sinπα∑_k=0^∞{(x/2)^α/Γ(k+1+α)[ log(x/2) - ψ(k+1+α) ]
    	+
    	
    (x/2)^-α/Γ(k+1-α)[ log(x/2) - ψ(k+1-α) ]
    },

we have ∂∂α K_α(x) = K̃_α(x) - π(πα) K_α(x), and so

    R_n (x) = x K_1-i ν_n(x)(x)/K̃_iν_n(x)(x).

Noting that from Lemma <ref> and 5.74 we have 

    ( x/2)^± i ν_1(x)   =
    	-1 ± i πγ/log(2/x) - γ + O( 1/|log x|^4), 
    	
    1/Γ(1 ± i ν_1)   = 
    	1 ± i γν_1(x) + O( 1/|log x|^2),
    	
    ψ(1 ± i ν_1 (x))
    	   =
    	-γ±ζ(2) ν_1(x) + O( 1/|log x|^2),

paying attention to cancellations, (<ref>) gives

    K̃_i ν_1(x)(x)
    	=
    	i L^2(x)/π + O( 1 ),

where the k = 1,2, … terms in the sum are bounded in the same way as the similar sum appearing in the proof of Lemma <ref>.
Similarly, noting that 

    1/Γ(-iν_1(x))  =  -iν_1(x) [1 - i γν_1(x) + O(ν_1^2(x)]  ,

(from 5.71) and expanding K_1-α(x) as in (<ref>), we have

    K_1 - i ν_1(x)(x)
    	=
    	1/x + O( 1/x |log x|^2),

taking advantage of a cancellation between the subleading terms in (x/2)^1+iν_1(x) and 1/Γ(1+iν_1(x)),
and so

    R_1(x) = - i π/[ log(2/x)-γ]^2 + O( 1/|log x|^4).





	For any a ∈ (1,2), There exist X_a, C_a > 0 such that 
	
    |
    		x K_1-α(x)/K_α(x)
    		+ 2 i R_1(x) ν_1(x)/α^2 + ν_1^2(x)|
    		< 
    		C_a/|log 2/x - γ|  ,

whenever x ∈ (0,X_a] and |α| = a (π / |log (2/x) - γ|).



From Lemma <ref> we see that we can choose 
X_a < 2 e^-γ such that 

    ν_2 (x)  >  a (π / |log (2/x) - γ|)  > ν_1(x)  ,

for all x ∈ (0,X_a]; then the quantity to be bounded is a continuous function of both α and x for all relevant values except x=0,
so we need only check that

    lim sup_x ↘ 0|
    	x L(x) K_1- α̃/ L(x)(x)/K_α̃/L( x)(x)
    	+ 2 i L(x) R_1(x) ν_1(x)/(α̃/L( x))^2 + ν_1^2(x)|  ,

is bounded uniformly for |α̃| = a π, 
where for brevity L(x) := log 2/x - γ; in fact we will show that both terms in the sum are suitably bounded.
In fact

    lim_x ↘ 0( x/2)^α̃/ log x  =   e^α̃  lim_x ↘ 0Γ( 1 + α̃/log x)
    	  =   1  ,

and with (<ref>) and (<ref>) this implies that

    lim_x ↘ 0K_α̃/L(x)(x) /L(x)  = e^α̃ - e^-α̃/2 α̃  .

Noting that limΓ(α̃/L( x)) / L( x) =1/α̃ we also have 

    x K_1-α̃/L(x)(x)
    	x ↘ 0∼
    	x L(x)/2 α̃(x/2)^α̃/L(x) - 1/Γ( α̃/L(x))∼ e^α̃;

then

    lim_x ↘ 0x L(x) K_1-α̃/ L(x)(x)/K_α̃/L(x)(x) = 2 α̃/e^2α̃ - 1,

which, recalling |α̃| = a π∈ (π,2π), is indeed uniformly bounded.

Recalling (<ref>), we have

    lim_x ↘ 0L^2 (x) R_1(x)  =   - i π  ,

and thus

    lim_x ↘ 02 i L(x) R_1(x) ν_1(x)/(α̃/L( x))^2 - ν_1^2(x) =  
    	2 π^2/α̃^2 - π^2  ,

which is also uniformly bounded in a suitable fashion. Hence (<ref>) is proven
and therefore also the proof of  Lemma <ref> is complete.


Noting that

    - 2 i R_1(x) ν_1(x)/α^2 + ν_1^2(x)
    	=
    	R_1(x)/α- i ν_1(x)
    	-
    	R_1(x)/α+ i ν_1(x),

the expression examined above is an analytic function of α in the interior of the circles under consideration apart from removable singularities,
and so using Lemma <ref> and (<ref>) we have

	For any a ∈ (1,2), There exist some C_a,X_a > 0 and a sequence of functions I_n:(0,X_a] → such that
	
    x K_1-α(x)/K_α(x)
    		=
    		∑_n=0^∞ I_n(x) α^n
    		-2 i R_1(x) ν_1(x)/α^2 + ν_1^2(x) ,

	whenever x ∈ (0,X_a] and |α| ≤ a (π / |log x|), and
	
    |I_n(x)| ≤ C (  |log (2/x) - γ|/a π)^n-1  ,

	for all n.

We then have 

    .
    	1/n!∂^n/∂α^n[ x K_1-α(x)/K_α(x)]
    	|_α=0
    	=
    	I_n(x) 
    	-
    	{( -1 )^n/2 2 i R_1(x)/( ν_1(x) )^n+1,     n  even
    
    		0,     n  odd..




From (<ref>) and (<ref>) we have that for any < 2 e^-1-γ there exist finite C_R, C_ν such that 

    | R_1(x) +  i π/L^2(x)|
    	≤C_R/L^4(x) and |
    	1/ν_1(x) - L(x)/π|
    	≤C_ν/π L^2(x)  ,

for all x ∈ (0,], and so 

    |
    	R_1(x)/( ν_1(x) )^n+1
    	+  i  L^n-1(x)/π^n|
    	= 
    	|
    	R_1(x) ( L(x)/π +
    	( 1/ν_1(x)- L(x)/π))^n+1
    	-  i  L^n-1(x)/π^n| 
    	
    ≤| R_1(x) ( L(x)/π)^n+1
    	-  i  L^n-1(x)/π^n|
    	+ | R_1(x) |1/π^n+1∑_m=1^n+1n+1m 
    	L(x)^n+1-3m C_ν^m
    	
    ≤C_R/π^n+1L^n-3(x)
    	+
    	(π + C_R/L^2(x))
    	1/π^n+1∑_m=1^n+1n+1m 
    	L(x)^n-3m-1 C_ν^m
    	.

We have 

    ∫_0^ L^n(x) dx
    	=
    	2 e^-γ∫_log(2/) - γ^∞ L^n e^-L dL
    	=
    	2 e^-γΓ(n+1, )  ,

for n ≥ 0,
where := log(2/) - γ for brevity (note > 0 since we have assumed < 2 e^-γ), and where Γ(n,) := ∫_^∞ t^n-1 e^-t dt is the upper incomplete Gamma function <cit.>.
For n<0, since we have assumed ≤ 2 e^-γ-1 we have L(x) ≥ 1 for x ∈ (0,), and thus

    0 ≤∫_0^ L^n(x) dx
    	≤  .

Note that by combining 8.82 and 8.101 we obtain

    Γ(n+1,)/Γ(n,)
    	=
    	n + ^n e^-/Γ(n,)≥
    	n + ≥ 
    	n  ,

which can be applied iteratively to obtain

    Γ(n+1,)/Γ(n+1-m,)≥n!/(n-m)!  .

for m ≤ n.  Then

    ∑_m=1^n
    		nm
    		C_ν^m    ∫_0^ L^n-3m-1(x)
    		
       ≤∑_m=1^(n-1)/3nm(n-3m)!/(n-1)!
    		C_ν^m
    		Γ(n,)
    		+
    		∑_m=(n-1)/3+1^n 
    		nm C_ν^m
    		
       =
    		n
    		∑_m=1^(n-1)/31/m!(n-3m)!/(n-m)!
    		C_ν^m
    		Γ(n,)
    		+
    		∑_m=(n-1)/3+1^n 
    		nm C_ν^m
    		
       ≤
    		n Γ(n,)
    		∑_m=1^∞(n-3m)!/(n-m)!C_ν^m/m!
    		+∑_m=0^n nm C_ν^m
    		
       =
    		27/2 (2n-3)(2n-6)( e^C_ν-1 ) Γ(n,)
    		+( C_ν +1 )^n   ,

using the observation that

    (n-m)!/(n-3m)!≥2/3 n ( 2n-3/3)( 2n - 6/3)  ,

for 1 ≤ m ≤ n/3.
Using this to bound the second term on the right-hand side of of Inequality (<ref>) and bounding the other two terms similarly,
we see that the integral in x from 0 to  of the right-hand side of Inequality (<ref>) admits a bound of order π^-n-1Γ(n,)/n^2 for large n.
We also have 

    | ∫_0^ I_n(x) dx |
    	≤∫_0^ |I_n(x)| d x 
    	≤2 C_a e^-γ/a^n-1Γ(n,)/π^n-1,

for any a ∈ (1,2),
and so the dominant behavior of the even Taylor coefficients ^(2n)(0)/(2n)! for n large is that of

    ( -1 )^n+12/π^2n∫_0^( log(2/x) - γ/π)^n-1 d x
    	    =  
    	4 e^-γ( -1 )^n+1Γ(2n,)/π^2n
        ∼  
    	4 e^-γ( -1 )^n+1(2n-1)!/π^2n  ,

noting Γ(n,) ∼Γ(n) = (n-1)! <cit.>,
while the symmetry noted in (<ref>) imposes that '(0) = 4 η and ^(2n+1)(0) = 0  for n =1,2, …. 

The proof of Theorem <ref> is therefore complete.





§ SCALING LIMIT OF MATRIX PRODUCT: PROOF OF THEOREM <REF>



As announced, we generalize  the set-up of  (<ref>)-(<ref>) in the sense that 
we prove 



Consider a family of positive random variables {Z^}_∈ (0, _0) such that
(Z^=y)=0 for every y and such that for some >0 and ∈ we have

    lim_↘ 0[Z^ -1 ]/  = 1/2 ^2 (1-)
         and    lim_↘ 0[(Z^ -1)^2 ]/  =  ^2   .
 
Assume moreover that for every c>0

    lim_↘ 01/( | Z^ -1|  >  c )   =  0 ,

and 

    lim sup_↘ 0|[1/Z^]-1/|  < ∞  .

Then if we consider the model (<ref>)-(<ref>) with the IID sequence {Z^ (n)}_n=1,2, …
generalized to an arbitrary IID sequence with common law satisfying (<ref>)-(<ref>), then  
(<ref>) and (<ref>) hold true. 




Theorem <ref> directly implies Theorem <ref>: the cases of two more classes of distributions are treated just before the proof. 
Note that with (<ref>) we are in reality just assuming the existence of the two limits and that the second limit is not zero.
The second assumption, i.e.  (<ref>),
 barely fails to be a consequence of 
(<ref>). 
The third assumption, i.e. (<ref>), is used to control the  amount of the   mass of Z^ that is  close to zero:
it is not difficult to realize that, given (<ref>), replacing  (<ref>) with 
the stronger condition

    lim_↘ 0[1/Z^]-1/  = 1/2( -1) ^2  ,
 
leads to very little loss of generality. Moreover, we have assumed that the law of Z^ ha no mass just to be sure that we do not fall into a pathological case for the theory of product of random matrices, but all we need is a condition 
that guarantees the existence of the limit in (<ref>)
and that the Markov chain associated to matrix product is ergodic: this is  true in greater generality    <cit.>.


Before giving the proof let us show two classes of examples to which Theorem <ref> applies:
 

  * The distribution chosen in <cit.> falls into the class

    _1^-y^-1_(0, _1)(y)   ,

with

    _1  = _1(,)   =    1+1/ +(1-)/^2+ o( 1/^2)  .

(→∞) is the parameter that tunes the strength of the disorder and Theorem <ref> can be applied 
by setting  =^-2: let us verify the hypotheses.
We compute for every ν
 
    [ ( Z_)^ν]  = _1(,)^ν/1+ ν/  = 1+ ν (ν-)/2 ^2 + o( 1/^2)  ,
 
and we directly obtain 

    lim_→∞^2[  Z_ -1
    ]   = 1-/2    and    lim_→∞^2[ ( Z_ -1)^2
    ]   =  1  ,

and 

    [ Z_^± 2]   =  
      1+  (2∓)/ ^2 + o( 1/^2)  = exp(   (2∓)/ ^2) 
      + o( 1/^2)   ,

Moreover or every c∈ (0,1) 
the event {|  Z_ -1| >c}= {  Z_ -1 <-c} if  is sufficiently large, because Z_≤_1(, ), which tends to one for →∞. On the other hand 
 ( Z_ -1 <-c)=((1-c)/ _1)^, which is bounded by (1-c)^ since _1>1.

  * Choose a centered and compactly supported probability density p(·) and set  ^2:=∫ t^2 p(t)  t. Then the  random variable Z^ with density given by 

    y ↦1/√() p ( y- m_/√())
         with  m_:= 1 + 1/2^2 (1-)   ,

with  smaller than a suitable _0>0, 
satisfies the hypotheses of Theorem <ref>.




Proof of Theorem <ref>.
We start with the proof of (<ref>), which  is a direct application of the approximation-diffusion principle: we exploit <cit.>, notably <cit.>.
Equivalently, one can resort to   <cit.>. The procedure demands three steps:


  * compute the local drift at ∈^2: uniformly for =(x_1, x_2)^ t in compact sets  

    b^()     = ^-1 A^ = 
    [           0            ;        [Z^] [Z^ -1 ]  / ]
       ↘ 0⟶   b()  :=  b  ,      with  b:= 
    [        0         ;          (1-)^2/2 ] ,

where we have applied the first assumption in (<ref>);

  * compute the diffusion matrix at : again uniformly we have

    a^()   = ^-1[A^ ^ t (A^)^ t]
    ↘ 0⟶  a() := [        0        0;        0 ^2 x_2^2 ]  ,

where we have applied both assumptions in (<ref>);

  * observe that, by  (<ref>), ^-1 (  |A^| ≥ c) → 0 for every c>0.




Then, since  the stochastic differential system  with drift b(·) and 
diffusion matrix a(·)   has  unique (strong) solution, the 
Markov chain X^ converges in law to the diffusion process
with drift b(·) and 
diffusion matrix a(·), which is precisely the solution X to the stochastic differential system
(<ref>).
This completes the proof of  (<ref>). 



In order to prove (<ref>) 
we start by observing that _Z^()=_Z^(-), in agreement with the analogous result for
_, (·)  (Theorem <ref>(2)), because D(I+A^)D, with
D the diagonal matrix with +1 and -1 on the diagonal, is equal to   I+A^ with  replaced by -.
Hence we can restrict to >0.
Moreover if
 we set
 Y^(n):= X_2^(n)/X_1^(n), we have  that Y^(⌊·/⌋) ⟶ Y(·) in law as ↘ 0 just because of  (<ref>) and because the map (x_1,x_2) ↦ x_2/x_1, from (0, ∞)^2 to 
 (0, ∞), is continuous. Denote by T_t^ and T_t the corresponding Markov operator semigroups

    T_t^ f(y)   =    E^_y[ f( Y^(⌊ t/⌋) ) ] ,     T_t f(y) =  E_y[ f( Y(t) ) ]  ,

acting on bounded continuous f: (0,∞) →. Note that we have also introduced the notation E^ and E for the expectation with respect
to the two Markov processes we consider.  We claim that:




  *  For bounded continuous f: (0,∞) → and t ∈ [0, ∞), we have that

    T_t^f(y) ↘ 0⟶ T_t f(y)    uniformly for  y  in  compact subsets of  (0,∞) .


  *  For all positive  there exists a unique law μ^ on (0,∞) which is  invariant  for the Markov chain Y^, which is ergodic. 

  *  
Choosing _0 ∈ (0, 1/) we have 

    sup_∈ (0,_0]
    ∫_0^∞ y^2μ^( y) < ∞ ,     
    sup_∈ (0,_0]∫_0^∞ y^-1μ^( y) < ∞ .




Claim (<ref>) is a byproduct of the proof of  (<ref>)  <cit.>.

Claim (<ref>) comes from the general theory of products of random matrices. Let us prove  (<ref>), and start by writing

    Y^(n+1) =   Z^(n+1)   u( Y^(n)) ,     u(y)= y+/1+ y .

Observing that for ∈ (0, 1/]

    ^2 / z^2(u( z^1/2 )^2) = -   (1-^2 ^2) 3z + 4  z^1/2  + 1/2 z^3/2 ( z^1/2 +1)^4≤ 0  ,

we obtain by the Markov property and by Jensen's inequality that for a given initial condition y>0

    E_y[ Y^(n+1)^2]    = [(Z^(n+1))^2)] E_y[u( Y^(n))^2] 
       ≤   q_,+^2 
     u( E_y[ Y^(n)^2]^1/2)^2  ,

where 

    q_,+  := √([(Z^)^2)])(<ref>)= 1+ (1-/2) ^2  + o(^2)
      .

Therefore if we set x_n:=E_y[ Y^(n)^2]^1/2 we have 
x_n+1≤ q_,+ u(x_n) which directly entails that x_n< ∞ for every n  and, since u(·) is bounded and concave increasing with u(0)>0,
 the application q_,+ u(·) has only one  positive fixed point that attracts every positive number. 
 The fixed point x^+_, () is easily computed:
 
    x^+_, ()  = 1/2( q_,+-1/+ √(( q_,+-1/)^2
     + 4 q_,+)) 
     
    ↘ 0∼1/2( 
      (1-/2)/ ^2
      + √(( (1-/2)/ ^2)^2
     + 4 
     ))
        .

 Therefore 
 lim sup_n x_n  ≤ x^+_, () 








and x^+_, () is bounded for ↘ 0.
Since {Y_n^}_n=0,1, … converges in law to the  random variable Y^_∞ that is distributed according to 
μ^, by standard measure theory argument we infer that
[ (Y^_∞)^2]=
∫_0^∞ y^2μ^( y) ≤ (x^+_, ())^2
 which 
 proves the first claim in (<ref>). 
 
 For the other claim in (<ref>) it is useful to note that Y^(n)=
Y^(n)^-1 evolves according to the similar dynamics driven by 1/Z^,

    Y^(n+1)  = (Z^(n+1))^-1Y^(n)+ /1+Y^(n)  .

We can now proceed in a simpler way than above and exploit directly the concavity of u(·) to get to

    x̃_n+1  :=  E_y[ Y^(n+1)] ≤   q_,- 
     u( E_y[ Y^(n)])  ,

and lim sup_n x̃_n ≤ x^-_,(), with x^-_,() defined replacing q_, + with 
q_, - in the definition (<ref>) of 
x^+_,(). It is therefore clear that
(<ref>) tells us that x^-_,() remains bounded for ↘ 0 and 
the second claim in (<ref>) is proven.



Of course if we make the stronger, but in practice almost equivalent, condition 
on the second moment of 1/Z^ in (<ref>), the argument for 
the first claim in (<ref>) applies and directly  yields 
sup_∈ (0,_0]∫_0^∞ y^-2μ^( y) < ∞.



With  (<ref>)–(<ref>) at hands, we complete the proof of  (<ref>).
By (<ref>)-(<ref>) and iterating we obtain

    log X_1^(n)     = log X_1^(n-1) + log(1+  Y^(n-1) )
      
       = log X_1^(0) + ∑_i=1^n log(1+  Y^(i-1) )    .

Following <cit.>,  we express the Lyapunov exponent

    _Z^ () = lim_n →∞1/nlogX^(n)  = lim_n →∞1/nlog X_1^(n)
    (<ref>)=lim_n →∞1/n∑_i=1^n log(1+  Y^(i-1) )   
    = ∫_0^∞log( 1+  y)  μ^ ( y)   .

By (<ref>), the family {μ^}_∈ (0,_0] of probability measures is tight on (0,∞).
By (<ref>) and <cit.>, every weak limit of {μ^}_∈ (0,_0] is invariant for Y, whose unique invariant measure has the density p_(·), implies that 
μ^( y) converges weakly to  p_(y)  y as ↘ 0. Then,

    _Z^ ()/ - _,  ()
     = ∫_0^∞(log(1+ y) / - y) μ^ (  y) +  ∫_0^∞  y  ( μ^ ( y)  -  p_(y) y)

The last term vanishes as ↘ 0 by weak convergence and  uniform integrability from claim (<ref>). But also  the first term in the right-hand side vanishes for the same reasons because

    |log( 1+  y) / - y |  = ∫_0^y ^2 z/1+ z z 
     ≤ √(∫_0^y  z
     ∫_0^y ^2 z z)  = ^3/2/√(2)^1/2 y^3/2 ,

where we have used that, for u≥ 0,  u/(1+u) is bounded above both by 1 and by u.
This completes the proof of (<ref>) and, therefore, also the proof of Theorem <ref>.






§ THE MCCOY-WU MODEL


In <cit.>, McCoy and Wu examined a two-dimensional Ising model with bond disorder of a particular type (subsequently known as the McCoy-Wu model): the couplings between sites in neighboring columns have a constant strength E_1, while the couplings between neighboring sites in the same column take  a random value E_2(n) which is fixed within each row but varies independently – keeping the same distribution – between different rows (Figure <ref>).  They showed that in the thermodynamic limit the free energy per site of this model is given (up to the subtraction of an analytic function of β) by  

    _MW()  := 1/4π∫_-π^π^MW_(θ) θ ,

where 
^MW_(θ) is the Lyapunov exponent of the random matrix

    M_β(θ)
    :=
    [         1 a/a^2+b^2; a/a^2+b^2  /a^2+b^2 ]  ,

with

    a(θ)  =  -2z_1 sin(θ)/| 1+z_1 exp(iθ) | ^2   and  b(θ)  = 1-z_1^2/| 1+z_1 exp(iθ) | ^2  ,
 
where 

    z_1  = tanh(  E_1)  ,   z_2(n)  = tanh(  E_2(n))
        and   =(n)= z_2^2(n)  .

In <cit.> a different version of the model has been considered: vertical bounds are random
in the horizontal direction and randomness is repeated in each line. This model, that allows frustration, is richer, but 
the features that are novel with respect to the McCoy-Wu model cannot be appreciated in the weak disorder limit: our analysis applies to <cit.> as well, but we will not develop this issue here. 




To avoid trivialities we assume that E_1≠ 0 as well as that E_2 is a non degenerate random variable:
it is immediately clear that the sign of E_2 does not matter and just a little thought reveals that the sign of E_1
is irrelevant too. Therefore we assume that  E_1∈ (0, ∞) and that E_2 is a random variable taking values in (0, ∞).
It is helpful (mostly to simplify the presentation) to assume that E_2 takes values in [E_2^-, E_2^+], with 0<E_2^- <E_2^+< ∞. 



Moreover one directly sees that a(·) is odd and b(·) is even, which yields
that ^MW_(·) is even: in fact  D_± M_(θ) D_±=
M_(-θ), with D_± the diagonal matrix with (+1,-1) on the diagonal. Therefore:

    _MW()  := 1/2π∫_0^π^MW_(θ) θ .


McCoy and Wu claim that for every υ∈ (0, π) – our focus is on υ small – the function

    ↦1/2π∫_υ^π^MW_(θ) θ .

is real analytic on (0, ∞). This can be proven by applying the main result in <cit.> (see also <cit.>). We sketch the argument here by considering separately the case
θ bounded away from 0 and π and the case of θ near π: with >0 small



  * For θ∈ [, π -] the matrix M_β(θ) (with positive entries) maps the closure of  the cone Q – here Q is first quadrant without the axes, that is the set of vectors with positive entries – to   Q∪{0}. More precisely, 
by the  hypothesis we have made on the suport of Z,  for every ∈ (0, π/2) and every ∈(0,1) the closure of Q is mapped  into a cone whose closure  is a subset of Q∪{0} and this subset is the same 
for every choice of θ∈ [, π -] and every ∈ [,1/]. This 
uniform cone property implies the real analyticity of 
↦^MW_(θ) with a convergence radius that is bounded away from zero uniformly in θ∈ [, π -] and ∈ [,1/].

  * For θ∈ [π -,π] we argue by observing first that

    a(π) = 0, 
    	  
    	b(π)
    	=
    	1-z_1^2/(1-z_1)^2
    	=
    	1+ tanhβ E_1/1 - tanhβ E_1
    	=
    	e^2 β E_1,

so

    M_β(π)
    	 = [                      1                      0;                      0 e^-4 β E_1tanh^2 β E_2 ].

Since e^-4 β E_1tanh^2 β E_2<1 the action of M_β(π) contracts uniformly any cone 
of the form {(x,y):   y ≥| x|}, in the sense there exists >0 such that  M_β(π)
sends {(x,y):   y ≥| x|} into {(x,y):   y ≥ (1+ ) | x|}, uniformly in >0 and E_2.
Elementary arguments show that this result is only slightly perturbed if we consider θ∈ [π -,π] with  sufficiently small. 
This uniform cone property implies the real analyticity of 
↦^MW_(θ) with a convergence radius that is bounded away from zero uniformly in θ∈ [π -,π] and >0.



Therefore the true issue is the regularity (or lack of it) of 

    ↦1/2π∫_0^υ^MW_(θ) θ ,

for a υ>0 that can be chosen as small as one wishes. 
At this point McCoy and Wu claim that 
the only non analytic point of the map in (<ref>) can be at _c defined by

    2 β_c E_1 +  [ logtanhβ_c E_2 ]   =   0  .


To see that this is the only possible candidate,
McCoy and Wu point out that

    a(0) = 0, 
    	  
    	b(0)  
    	=
    	1-z_1^2/(1+z_1)^2
    	=
    	1- tanhβ E_1/1 + tanhβ E_2
    	=
    	e^-2 β E_1,

so 

    M_β(0)
    	=
    	[                     1                     0;                     0 e^4 β E_1tanh^2 β E_2 ],

and so

    ^MW_β(0)
    	=
    	max( 0, 4 β E_1 + 2 [ logtanhβ E_2 ] )

for β real. This admits an analytic extension in a neighborhood of any positive β except for β_c.

This is of course far from being close to a proof, since one has to control the integral over θ∈ (0, υ) and
not the value in zero. But McCoy and Wu  perform also a more subtle analysis that can be understood 
precisely via the diffusion limit of matrix products that is at the center of our analysis. To explain this  let us make a further manipulation 
to match more sharply our framework.



In fact, as it stands, 
M_β(θ), cf. (<ref>), is not of the form 
(<ref>). But by 
noting that 

    1/a^2(θ)+b^2(θ)
    	=
    	(1+z_1)^4/(1-z_1^2)^2 + O(θ^2)
    	=
    	( 1+z_1/1-z_1)^2+ O(θ^2)
    	=
    	e^4 β E_1+ O(θ^2)  ,

and

    a/a^2(θ)+b^2(θ) 
    	= 
    	-2 ( 1+z_1/1-z_1)^2 θ + O(θ^2) = -2
    	e^4 β E_1+ O(θ^2)  ,


if we let

    ε̃:= 
    	2 z_1/(1-z_1)^2θ
    	
    	
    	
    	
    	  ,


we see  that to leading order as θ↘ 0

    [         1       -ε̃;       -ε̃ e^4 β E_1 ]

is  M_β(θ). The matrix in (<ref>)
is of the form (<ref>) up to a conjugation and a change of variables: 
in fact 

    [ 1  ; Z Z ] := [          1 e^-2 β E_1; e^-2 β E_1  e^4 β E_1 ]  = [           1           0;           0 - e^2 β E_1 ][         1       -ε̃;      - ε̃ e^4 β E_1 ][            1            0;            0 - e^-2 β E_1 ]  ,

and we observe – recall (<ref>) – that
 = c_θ, with c_= 2 sinh(2 E_1).




It is important to remark at this stage that the inverse temperature  and our fundamental parameter
 – we recall that  is the unique non zero real solution to
Z^=1 (Z=e^4 E_1tanh^2( E_2) depends on !) when such a solution exists and otherwise =0 – should be seen as an analytic change of variable: this is treated in detail in Lemma <ref>. In particular 
(_c)=0 and therefore
()= (-_c)  ' (_c)+ O((-_c)^2), but the constant ' (_c) depends of the law of Z (with =_c) and this expansion
should be done more carefully when the disorder is weak because, as we will see, ' (_c) becomes large
in this limit: this is treated in (<ref>)-(<ref>).



What McCoy and Wu do at this point is 


  * making a specific choice of 
Z=Z^=Z^_ that satisfies the hypotheses of Theorem <ref> (say, with =1 for simplicity); this  actually implements two choices:

 
  * the first is evident and it is the fact that disorder can be made weak by making  small;
 
  * the second is that -_c is chosen small and, precisely, of the order of . As we will
 explain, if we set y=(-_c)/ and we keep y∈ fixed, then () ∼ -C__c y, and the constant C__c>0
will be given explicit in the specific case that we are going to develop, see (<ref>).


  * they choose also υ∝: let us fix in an arbitrary fashion υ =.



In physical terms these choices correspond to focusing on the critical window in the limit of weak disorder. Cutting
the integral at θ= is harmless (as we have discussed before), but of course only as far as  is kept fixed.

McCoy and Wu  are in the end just dealing (recall (<ref>)) with the Lyapunov exponent _, _c + y (c__cx )
of the matrix (we perform the change of variable 
θ=x)

    [                 1            c__c x; c__c x  Z__c + y^         Z__c + y^ ]  .
 
But Theorem <ref> (see also Theorem <ref>) tells
us that _, _c + y (c__c x) is asymptotically equivalent for  small to _1, C__c( c__cx) so that

    ∫_0^^MW__c + y (θ) θ ∼ ∫_0^1 _, _c + y (c__c x) x ∼ ^2 ∫_0^1 _1,C__c( c__cx) x

and we remind the reader that _1,C__c( c__cx) has the explict expression (<ref>).
Therefore, up to two inessential  constants we arrived at (<ref>). We did not fully justify the 
equivalences in (<ref>), but this is not really the main problem: the main unresolved mathematical issue is that 
what we are after is proving that, for a fixed (possibly extremely small) value of ,  the leftmost term in 
(<ref>) is a C^∞ function of y at 0 and that the same expression is not analytic at zero. McCoy and Wu instead argue (and we prove in Theorem <ref>) that 
↦∫_0^1 _1,C__c( c__cx) x has these properties: but this second statement does not imply the first. 


We now complement our discussion with the analysis of the specific distribution chosen for the disorder law in <cit.>. We also 
discuss more in detail the change of variable (). 



  §.§.§ Analysis of the distribution chosen by McCoy and Wu <cit.>

 McCoy and Wu consider the disordered variable =tanh^2(  E_2) that depends on a parameter that they call  and it is large: in fact
 
    = ^-2  .

 The density  of  is 
  supported on (0, _0)   and  equal to 
 _0^-y^-1 for y∈ (0, _0). Necessarily  _0=_0()= tanh^2( E_2^*), with 
 E_2^* the maximum value that the random variable E_2 can reach. The density of Z=  Z_ 
 (recall that Z is defined in (<ref>))
 is therefore
 
    _1^-y^-1_(0, _1)(y)     with  _1()  =  e^4 E_1_0()  .

 Note that for every ν∈ (-, ∞) 
 
    [ ( Z_)^ν]  = _1()^ν/1+ ν/  ,
 
 and we want to solve for =()≠ 0 the equation
 
    [ (  Z_)^]  -1  = _1()^ - 1- //1+ /  =  0  .
 
  On one hand  we compute
 
    log (_1())  =
    logtanh ^2 (_c E_2^*) - [logtanh ^2 (_c E_2)] + 4 ( -_c)+ 
     logtanh ^2 ( E_2^*)- logtanh ^2 (_c E_2^*)  ,
 
 and a straightforward computation yields 
 
    [logtanh ^2 ( E_2)]  = logtanh ^2 ( E_2^*) - 1/  ,

 so  for  close _c we have 
 
    log_1()  = 1/+ ( -_c) ( 4E_1 + 1/sinh (2 _c E_2^*)) + O( ( - _c)^2)  .
 
 On the other hand from (<ref>)
we see that if  is fixed (so we look at  as a function of ) we have 
 
    log_1()  = log(1+ /)/→∞= 1/- /2^2+ O( 1/^3)  .

By comparing (<ref>) and (<ref>) we see that if (-_c) ^2=O(1) then 

    ()   =  -( -_c) ^2( 8E_1 + 2/sinh (2 _c E_2^*)) + O( ^-1)  .




  §.§.§ On the relation between  and 

Here are the details of the important map that relates  and : 




Assume that the support of the random variable E_2 is bounded away from zero, so Z= exp(4  E_2) 
tanh^2( E_2) is supported on a compact subinterval of (0, ∞). Assume also that E_2 is not constant. 
Then the  equation

    [ Z^]-1/  =  0  ,

has a unique real solution  for every >0. This defines a map
↦() from (0, ∞) to . This map is decreasing, hence it is a bijection,  and it is real analytic.




Proof. 
Let f: × (0,∞) → be the function defined by
f(, ):=[ Z^]-1/, for ∈∖{0}
for α≠ 0, and f(0, ):=  [log Z].
It is straightforward to see, using the support properties of E_2, that f is real analytic on its entire domain. 
Then we observe that, for fixed , Z^ is an increasing function of  and, by the support properties,
this implies that ∂_ f(, )>0 for every >0 and ∈.
On the other hand  
 if we set g_ ()= [ Z^]-1
we have that ∂_ f(, )=(  g'_()- g_ ()) / ^2. 
But g_ (·) is  (strictly) convex and g_ (0)=0: so g'_()- g_ ()>0 for ≠ 0
and therefore ∂_ f(, )>0 for ≠ 0. For =0 it suffices to perform a Taylor expansion
of g_() at =0 to see that ∂_ f(, )|_=0= g^''_ (0)/2>0.
The proof is completed by  applying the  Implicit Function Theorem for real analytic functions <cit.>.





§ ACKNOWLEDGMENTS
 
We are very grateful to  Bernard Derrida for very insightful discussions, and to an anonymous referee for pointing out significant errors in the statement and proof of Proposition <ref> in an earlier manuscript. 
Part of this work was developed while G. G. and R. L. G. were visiting IHP (Paris) during the spring-summer 2017 trimester.  G. G. acknowledges the support of grant ANR-15-CE40-0020.  The work of
R. L. G.was funded by the European Research Council under the European Union's Horizon 2020 Programme, ERC Consolidator Grant UniCoSM (grant agreement no 724939).


99








cf:AKQ
Alberts, T.,  Khanin, K. and Quastel, J.; The intermediate disorder regime for directed polymers in dimension 1 + 1. Ann. Probab. 42 (2014), 1212-1256.

cf:Baxter 
 Baxter, R. J.; Exactly solved models in statistical mechanics. Academic Press,  1982.

cf:BL
   Bougerol, P.  and Lacroix, J.;
    Products of random matrices with applications to Schrödinger
              operators.
   Progress in Probability and Statistics 
    8,
 Birkhäuser Boston, Inc., Boston, MA,
      1985.
      
cf:bouchaud 
 Bouchaud, J.-P., Comtet, A., Georges, A.  and Le Doussal, P.;
Classical diffusion of a particle in a one-dimensional random force field, 
Ann. Physics 201 (1990),   285-341.       
     
     
cf:CSZ      
 Caravenna, F.,   Sun, R. and  Zygouras, N.;
Polynomial chaos and scaling limits of disordered systems.
J. Eur. Math. Soc. 19 (2017), 1-65.

CattiauxChafaiGuillin
Cattiaux, P.,  Chafaï, D. and Guillin, A.; 
Central limit theorems for additive functionals of ergodic Markov diffusions processes.
ALEA Lat. Am. J. Probab. Math. Stat.  9 (2012), 337-382.

      
  cf:contemp
   Cohen, J. E.,   Kesten, H.  and  Newman,  C. M. (Editors); 
   Random matrices and their applications. 
Proceedings of the AMS-IMS-SIAM joint summer research conference held at Bowdoin College, Brunswick, Maine, June 17-23, 1984. Contemporary Mathematics, 50, American Mathematical Society,  1986.

cf:desbois
 Comtet, A., Desbois, J. and  Monthus, C.;  Localization properties in one-dimensional disordered supersymmetric quantum mechanics, Ann. Physics 239 (1995),  312-350. 


cf:CLTT
Comtet, A., Luck, J. M., Texier, C. and
 Tourigny, Y.;
 The Lyapunov Exponent of Products of Random 2× 2 Matrices Close to the Identity.
  J. Statist. Phys. 150 (2013), 13-65.
  
 cf:CTT Comtet, A., Texier, C. and  Tourigny, Y.;
Lyapunov exponents, one-dimensional Anderson localization and products of random matrices, 
J. Phys. A 46 (2013), 254003, 20 pp.  


CPV
Crisanti, A., Paladin, G. and Vulpiani, A.;
Products of random matrices in statistical physics, volume
  104 of Springer Series in Solid-State Sciences.
Berlin: Springer-Verlag, 1993.


 
 
DH
Derrida, B. and Hilhorst, H. J.;
Singular behaviour of certain infinite products of random 2× 2 matrices.
J. Phys. A 16 (1983), 2641-2654. 


  
  cf:dubois
   Dubois, L.;
    Real cone contractions and analyticity properties of the characteristic exponents,  Nonlinearity 21 (2008),  2519-2536.

Dunster
Dunster, T.M.;
Bessel functions of purely imaginary order, with an application to second-order linear differential equations having a large parameter.
SIAM J. Math. Anal. 21 (1990), 995-1018.

EthierKurtz
Ethier, S. and Kurtz, T.;
Markov processes. Characterization and convergence.
Wiley Series in Probability and Mathematical Statistics: Probability and
   Mathematical Statistics.
1986.

cf:dfisher Fisher, D. S.;
 Critical behavior of random transverse-field Ising spin chains, Phys. Rev. B 51 (1995), 6411-6461.
 

cf:FW
 Freidlin, M. I. and Wentzell, A. D.;  Random perturbations of dynamical systems.  Third edition. Grundlehren der Mathematischen Wissenschaften, 260,  Springer, 2012.

Friedlander
  Friedlander, F. G.;
  Diffraction of Pulses by a Circular Cylinder,
  Comm. Pure and Appl. Math. 7 (1954), 705-732.
  
 cf:FrLl
 Frisch, H. L. and Lloyd, S. P.;
 Electron levels in a one-dimensional random lattice, Phys. Rev. 120 (1960), 1175-1189.
 

cf:GGG
Genovese, G.,  Giacomin, G. and  Greenblatt, R. L.;
Singular behavior of the leading Lyapunov exponent of a product of random 2×2 matrices,
Commun. Math. Phys. 351 (2017), 923-958. 

cf:GB Giacomin, G.; Random polymer models, Imperial College Press, World Scientific, 2007.
 
 cf:G Giacomin, G.;
  Disorder and critical phenomena through basic probability models, École d'été de probablités de Saint-Flour XL-2010, Lecture Notes in Mathematics 2025, Springer, 2011.

cf:grabsch 
Grabsch, A., Texier, C. and   Tourigny, Y.; One-dimensional disordered quantum mechanics and Sinai diffusion with random absorbers,  J. Stat. Phys. 155 (2014),  237-276.

cf:benjamin B. Havret, Regular expansion for the characteristic exponent of a product of 2 × 2 random matrices, arXiv:1804.06166

cf:KS 
Karatzas, I. and  Shreve, S. E.;
Brownian Motion and Stochastic Calculus.
Graduate Texts in Mathematics
113, Springer, 1998.

cf:primer
 Krantz, S. G. and  Parks, H. R.;
 A primer of real analytic functions,   Birkhüser Verlag,  1992.

Kaup
Kaup, L. and Kaup, B.;
Holomorphic Functions of Several Variables.
1983.






cf:Luck
 Luck, J. M.;   Critical behavior of the aperiodic quantum Ising chain in a transverse magnetic field, J. Statist. Phys. 72 (1993),  417-458.

MaruyamaTanaka
Maruyama, G. and Tanaka, H.;
 Ergodic property of N-dimensional recurrent Markov processes. Mem. Fac. Sci. Kyushu Univ. Ser. A  13 (1959), 157-172.

MWbook
McCoy, B. M. and Wu, T. T.;
The Two-Dimensional Ising Model. Harvard University Press, 1973.

MW1
McCoy, B. M. and Wu, T. T.;
Theory of a two-dimensional Ising model with random impurities.
  I. Thermodynamics.
Phys. Rev. 176 (1968), 631-643.

DLMF
Olver, F. W. J., Olde Daalhuis, A. B.,  Lozier, D. W.,  Schneider, B. I.,  Boisvert, R. F.,   Clark, C. W.,   Miller, B. R. and  Saunders, B. V. (Editors);
  NIST Digital Library of Mathematical Functions.
  http://dlmf.nist.gov/, Release 1.0.16 of 2017-09-18. 


NL
 Nieuwenhuizen, Th. M. and Luck, J. M.;  Exactly soluble random field Ising models in one dimension. 
J. Phys. A 19 (1986), 1207-1227.

PardouxVeretennikov
 Pardoux, E. and  Veretennikov, A. Yu.;
On the Poisson equation and diffusion approximation. I. Ann. Probab. 29 (2001), 1061-1085.

cf:ramola
 Ramola, K. and Texier, C.; Fluctuations of random matrix products and 1D Dirac equation with random mass, J. Stat. Phys. 157 (2014), 497-514.

RogersWilliams
Rogers, L. and Williams, D.;
Diffusions, Markov processes, and martingales. Vol. 2. 
Itô calculus. Cambridge University Press, Cambridge, 2000.

Ruelle
Ruelle, D.;
Analytic properties of the characteristic exponents of random matrix products.
Adv. Math 32 (1979), 68-80.

cf:Ruelle-book
Ruelle, D.; Thermodynamic formalism. The mathematical structures of equilibrium statistical mechanics. Second edition. Cambridge Mathematical Library. Cambridge University Press,  2004.

cf:Sadel 
Sadel, C. and  Schulz-Baldes, H.; Random Lie group actions on compact manifolds: a perturbative analysis.
Ann. Probab. 36 (2010), 2224-2257.

cf:schomerus
 Schomerus, H.,  and  Titov, M.;
Statistics of finite-time Lyapunov exponents in a random time-dependent potential, 
Phts. Rev. E 66 (2002), 066207 (11 pages).

ShankarMurthy
Shankar, R. and Murthy, G.;
Nearest-neighbor frustrated random-bond model in d=2: Some
  exact results.
Phys. Rev. B 36 (1987), 536-545.


StroockVaradhan
Stroock, D. and Varadhan, S.; 
Multidimensional diffusion processes.
Classics in Mathematics. Springer-Verlag, Berlin, 2006.

cf:marina Taleb, M.; 
Large deviations for a Brownian motion in a drifted Brownian potential, 
Ann. Probab. 29 (2001), 1173-1204. 

cf:Zanon
Zanon, N. and  Derrida, B.; 
Weak disorder expansion of Lyapunov exponents in a degenerate case. J. Stat.
Phys. 50 (1988), 509-528.



