Nonlinear $q$-Stokes phenomena for $q$-Painlev\'{e} I

We consider the asymptotic behaviour of solutions of the first $q$-difference Painlev\'{e} equation in the limits $|q|\rightarrow 1$ and $n\rightarrow\infty$. Using asymptotic power series, we describe four families of solutions that contain free parameters hidden beyond-all-orders. These asymptotic solutions exhibit Stokes phenomena, which is typically invisible to classical power series methods. In order to investigate such phenomena we apply exponential asymptotic techniques to obtain mathematical descriptions of the rapid switching behaviour associated with Stokes curves. Through this analysis, we also determine the regions of the complex plane in which the asymptotic behaviour is described by a power series expression, and find that the Stokes curves are described by curves known as $q$-spirals.


Introduction
In this paper we study the first q-difference Painlevé equation (q-P I ) where w = w(x), w = w(qx), w = w(x/q), and q ∈ C such that |q| = 0, 1. In particular, we assume that |q| > 1 and consider (1) under the limits |q| → 1 and n → ∞. Equation (1) is part of a class of integrable, second-order nonlinear difference equations known as the discrete Painlevé equations that tend to the ordinary Painlevé equations in the continuum limit. More generally, equation (1) is known as a q-difference equation since the evolution of the independent variable x takes the form x = x 0 q n for some initial x 0 . Motivated by Boutroux's study of the first Painlevé equation [7], the asymptotic behaviour of (1) in the limit |x| → ∞ has been considered in [28,30]. Joshi [28] showed that there exists a true solution satisfying w → 0 as |x| → ∞, which is asymptotic to a divergent series. However, [28] does not describe the Stokes switching behaviour that is typically associated with (divergent) asymptotic series.
In this study we uncover the Stokes behaviour present in the solutions of (1) using exponential asymptotic methods. We first introduce a parameter, , such that the double limit |q| → 1 and n → ∞ is equivalent to → 0. Since we may parametrize x by x = x 0 q n , one possible scaling choice which captures the desired behaviour is n = s/ and q = 1 + . Under this choice of scaling we have x ∼ x 0 e s + O( ) as → 0. We note that the large x limit may be obtained if s is taken to be sufficiently large and positive.
Exponential asymptotic techniques for differential-difference equations were developed by King and Chapman [35] in order to study a nonlinear model of atomic lattices based on the works of [13,48]. The authors of [31,32] applied the Stokes smoothing technique described in [35] to the first and second discrete Painlevé equations and obtained asymptotic approximations which contain exponentially small contributions. The difference equations considered in [31,32] are of additive type as the independent variable is of the form z n = αn + β. Motivated by their work, we extend this to q-P I in order to study asymptotic solutions of q-difference equations which display Stokes phenomena.
We note that there are other exponential asymptotic approaches used to study difference equations [25,47,49]. In particular, Olde Daalhuis [47] considered a certain class of second-order linear difference equations, and applied Borel summation techniques in order to obtain asymptotic expansions with exponentially small error. The methods used in these studies involves applying Borel resummation to the divergent asymptotic series of the problem, which allows the optimallytruncated error to be computed directly. We note that Borel summation techniques may be applied whenever a series possess factorial-over-power divergence in the late-order terms and therefore could be applied as an alternative approach to that utilised in the present study.

Background
General solutions to the Painlevé equations are higher transcendental functions, which cannot be expressed exactly in terms of elementary functions and therefore much of the analytic information regarding their general solutions are very limited. In fact, Nishioka [45] proved that the general solutions of q-P I are not expressible in terms of solutions of first-order q-difference equations.
The discrete Painlevé equations have appeared in areas of physical interest such as in mathematical physics models describing two-dimensional quantum gravity [10,20,22,51,52,59]. The basis of these models have origins in orthogonal polynomial theory, in which discrete Painlevé equations have been found to commonly appear [36,[38][39][40][55][56][57]. Since the discrete Painlevé equations often arise in many nonlinear models of mathematical physics they are often regarded as defining new nonlinear special functions [15,26].
Motivated by these applications, much research has gone into the asymptotic study of the (discrete) Painlevé equations. Previous asymptotics studies for the first discrete Painlevé equation have been conducted in [27,31,58] where the authors found solutions asymptotically free of poles in the large independent variable limit, which share features with the (tri)-tronqueé solutions of the first Painlevé equation found by Boutroux [7]. Asymptotic solutions have also been found for the so called alternate discrete Painlevé I equation [16,34].
Extending the work of [31], the authors of [32] also find solutions of the second discrete Painlevé equation (dP II ), which are asymptotically free of poles in some domain of the complex plane containing the positive real axis. This was achieved by rescaling the problem such that the step size of the rescaled difference equation is small in the asymptotic limit. A similar study on dP II was also investigated by Shimomura [54]. Shimomura showed that an asymptotic solution of dP II which reduces to the tri-tronqueé solution of P II can be found by choosing a scaling of dP II such that the rescaled equation tends to the second continuous Painlevé equation (P II ) in the limit n → ∞.
The isomonodromic deformation method has also been used to study the asymptotics for nonlinear difference equations [11,21,37,59]. Using the nonlinear steepest descent method developed by Deift and Zhou [17], the authors of [60] find the asymptotic behaviour of solutions of the fifth discrete Painlevé equation in terms of the solutions of the fifth continuous Painlevé equation.
However, there have been very few asymptotic studies on q-Painlevé equations. The asymptotic study of variants of the sixth q-Painlevé equation have been investigated by Mano [41] and Joshi and Roffelson [33] using the q-analogue of the isomonodromic deformation approach. The first q-Painlevé equation has also been investigated in [28,30]. In particular, Joshi [28] proves the existence of true solutions of (1), which are asymptotic to a divergent asymptotic power series in the limit |x| → ∞. Such expansions are known to exhibit Stokes phenomena in the complex plane.
Stokes phenomena are generally well understood in the case of linear and nonlinear differential equations [9,12,14,24,50]. The study of Stokes phenomena have also been investigated for nonlinear difference equations [35], and have been extended to study discrete Painlevé equations [31,32].
In the continuous theory, Borel summation methods are also used to describe Stokes behaviour by resumming divergent asymptotic series expansions. The q-analogues of these methods have also been developed for linear q-difference equations [18,53,[61][62][63] in order to describe behaviour known as q-Stokes phenomena. These methods have been explicitly applied to certain classes of second-order linear q-difference equations by Morita [42][43][44] and Ohyama [46].
However, q-Stokes phenomenon is unlike classical Stokes phenomenon for differential equations. The notion of q-Stokes phenomenon and its differences to classical Stokes phenomenon are detailed in [18,53]. To the best of our knowledge there have been no corresponding studies for nonlinear q-difference equations. The goal of this study is to extend the exponential asymptotic methods used in [31,32,35] to q-difference equations in order to describe Stokes behaviour present in the asymptotic series expansions.

Exponential asymptotics and Stokes curves
Conventional asymptotic power series methods fail to capture the presence of exponentially small terms, and therefore these terms are often described as lying beyond-all-orders. In order to investigate such terms, exponential asymptotic methods are used. The underlying principle of these methods is that divergent asymptotic series may be truncated so that the divergent tail, also known as the remainder term, is exponentially small in the asymptotic limit [8]. This is known as an optimally-truncated asymptotic series. Thereafter, the problem can be rescaled in order to directly study the behaviour of these exponentially small remainder terms. This idea was introduced by Berry [3][4][5], and Berry and Howls [6], who used these methods to determine the behaviour of special functions such as the Airy function.
The basis of this study uses techniques of exponential asymptotics developed by Olde Daalhuis et al. [48] for linear differential equations, extended by Chapman et al. [13] for application to nonlinear differential equations, and further developed by King and Chapman [35] for nonlinear differential-difference equations. A brief outline of the key steps of the process will be provided here, however more detailed explanation of the methodology may be found in these studies.
In order to optimally truncate an asymptotic series, the general form of the coefficients of the asymptotic series is needed. However, in many cases this is an algebraically intractable problem. Dingle [19] investigated singular perturbation problems and noted that the calculation of successive terms of the asymptotic series involves repeated differentiation of the earlier terms. Hence, the late-order terms, a m , of the asymptotic series typically diverge as the ratio between a factorial and an increasing power of some function as m → ∞. A typical form describing this is given by the expression as m → ∞ where Γ is the gamma function defined in [1], while A, and χ are functions of the independent variable which do not depend on m, known as the prefactor and singulant respectively. The singulant is subject to the condition that it vanishes at the singular points of the leading order behaviour, ensuring that the singularity is present in all higher-order terms. Chapman et al. [13] noted this behaviour in their investigations and utilize (2) as an ansatz for the late-order terms, which may then be used to optimally truncate the asymptotic expansion.
Following [48] we substitute the optimally-truncated series back into the governing equation and study the exponentially small remainder term. When investigating these terms we will discover two important curves known as Stokes and anti-Stokes curves [2]. Stokes curves are curves on which the switching exponential is maximally subdominant compared to the leading order behaviour. As Stokes curves are crossed, the exponentially small behaviour experiences a smooth, rapid change in value in the neighbourhood of the curve; this is known as Stokes switching. Anti-Stokes curves are curves along which the exponential term switches from being exponentially small to exponentially large (and vice versa). We will use these definitions to determine the locations of the Stokes and anti-Stokes curves in this study.
By studying the switching behaviour of the exponentially small remainder term in the neighbourhood of Stokes curves, it is possible to obtain an expression for the remainder term. The behaviour of the remainder associated with the late-order terms in (2) typically takes the form SA exp(−χ/ ), where S is a Stokes multiplier that is constant away from Stokes curves, but varies rapidly between constant values as Stokes curves are crossed. From this form, it can be shown that Stokes lines follow curves along which χ is real and positive, while anti-Stokes lines follow curves along which χ is imaginary. A more detailed discussion of the behaviour of Stokes curves is given in [2].

Paper outline
In Section 2, we find two classes of solution behaviour of q-P I , which we refer as Type A and Type B solutions. We also determine the formal series expansions of these solutions, and provide the recurrence relations for the coefficients.
In Section 3, we determine the form of the late-order terms for Type A solutions and use this to determine the Stokes structure of these asymptotic series expansions. We then calculate the behaviours of the exponentially small contributions present in these solutions as Stokes curves are crossed. This is then used to determine the regions in which the asymptotic power series are accurate representations of the dominant asymptotic behaviour.
In Section 4, we consider Type B solutions of q-P I following the analysis in Section 3. In Section 5 we establish a connection between Type A and B solutions to the vanishing and non-vanishing asymptotic solutions of q-P I found by Joshi [28]. Finally, we discuss the results and conclusions of the paper in Section 6. Appendices A-D contain detailed calculations needed in Section 3.

Asymptotic series expansions
In this section, we expand the solution as a formal power series in the limit → 0, obtain the recurrence relation for the coefficients of the series and deduce the general expression of the lateorder terms.
We first rewrite (1) as an additive difference equation by setting x = x 0 q n , which gives where w n = w(x 0 q n ). In our analysis, we will introduce a small parameter, , by rescaling the variables appearing in (3). The choice of scalings we apply are given by Under these scalings, equation (3) becomes and consider the limit → 0. We also note that under the scalings given by (4), the independent variable, x, has behaviour described by as → 0. As e s is an entire function, we set x 0 = 1 for the remainder of this analysis. It can be shown that equation (5) is invariant under the mapping as → 0, and where λ 3 = 1. We note that this corresponds to the rotational symmetry of (1) found in [28]. We expand the solution, W (s), as an asymptotic power series in by writing as → 0. Substituting (8) into (5) and matching terms of O( r ) we obtain the nonlinear recurrence relation for r ≥ 0 and where the polynomials P n (s) are given by where s 1 (n, k) are the Stirling numbers of the first kind [1]. This recurrence relations allows us to calculate W r in terms of the previous coefficients. From (9) we find that the leading order behaviour satisfies Equation (10) is invariant under s → s + 2πi as the function e s is 2πi-periodic and hence W 0 (s) is 2πi-periodic. Furthermore, it will be shown in Section 3.2 that the Stokes switching behaviour of (8) depends on W 0 (s) to leading order. Hence we restrict our attention the domain, D 0 , described by We also define the domain, D k , which we call k th -adjacent domain, by As W 0 satisfies a quartic we therefore have four possible leading order behaviours as → 0. We first define the following and Then the four solutions for W 0 are given by and From (15) and (16) it can be shown that Each of the leading order solutions, W 0 , are singular at points for which the argument of the square root term of B is equal to zero. The singularities of W 0 are given by where k ∈ Z. Let us denote the singularities in D 0 by Then it can be shown that the local behaviour of W 0,j (s) near the singular points (18) is given by as s → s 0,j for j = 1, 2, 3. Equation (17) shows that W 0,4 (s) is the sum of W 0,j (s) for j = 1, 2, 3, and is therefore singular at the points s 0,1 , s 0,2 and s 0,3 in D 0 . Two types of leading order behaviours can be characterized by the number of points at which they are singular in D 0 . In the subsequent analysis, we will refer to solutions with leading order behaviour described by W 0,j for j = 1, 2, 3 as Type A, while those with leading order behaviour described W 0,4 as Type B.
In Sections 3.3 and 4.1 we will show that the (anti-) Stokes curves emanate from these singularities. Consequently, the Stokes structure of Type A solutions will be shown to emerge from a single singularity in D 0 , while Type B solutions will have (anti-) Stokes curves emanating from three singular points in D 0 . In this sense, Type B solutions will have the most complicated Stokes behaviour due to possible interaction effects between the distinct singularities.

Type A Exponential Asymptotics
In this section we will investigate the Stokes phenomena exhibited in Type A solutions with leading order behaviour W 0,3 as → 0. This is done by first determining the leading order behaviour of the late-order terms, which will allow us to optimally truncate (8) and study the optimally-truncated error. The results for the remaining Type A solutions may be obtained using the symmetry (7).

Late-order terms
As discussed in Section 1.2, the ansatz for the late-order terms is given by a factorial-over-power form since the determination of W r in (9) involves repeated differentiation. We therefore assume that the coefficients of (8) with W 0 = W 0,3 are described by as r → ∞, where χ 3 (s) is the singulant, U 3 (s) is the prefactor and γ 1 a constant. We substitute (21) into (9) to obtain where the remaining terms are negligible for the purposes of this demonstration as r → ∞. By matching terms of O(W r ) as r → ∞, we find that the leading order equation is given by as r → ∞. Matching at the next subsequent order involves matching terms of O(W r−1 ) in (22). Doing so gives as r → ∞. In order to determine the singulant, χ 3 (s), we consider (23). The leading order behaviour of χ 3 (s) may be determined by replacing the upper limit of the sum in (23) by infinity. Taking the series to be infinite introduces error in the singulant behaviour which is exponentially small in the limit r → ∞, which is negligible here [31,32,35]. Evaluating the sum appearing in (23) as an infinite series gives The solution of (25) is given by where M ∈ Z. Noting that there are two different expressions for the singulant, we name them χ 3 (s; M ) and χ − 3 (s; M ) with the choice of the positive and negative signs respectively. In general, the behaviour of W r will be the sum of expressions (21), with each value of M and sign of the singulant [19]. However, this sum will be dominated by the two terms associated with M = 0 as this is the value for which |χ| is smallest [13]. Thus, we consider the M = 0 case in the subsequent analysis. Hence we set In order to determine the prefactor associated with each singulant we solve (24). As before, we replace the upper limit of the sum appearing in (24) by infinity and use (26) in order to obtain It can be verified that where Υ is a constant of integration, solves the differential equation (28) without the 3W 1 /W 0,3 term. By setting U 3 (s) = F (s)φ(s) and substituting this into (28) we obtain the differential equation where we have used (25) in order to express (30) in terms of χ 3 . Integrating (30) then gives whereΥ is a constant of integration. Hence, if we let U 3 and U − 3 denote the prefactors associated with the singulants χ 3 and χ − 3 , respectively, then the solutions of (28) are given by where Λ,Λ are constants (which depend on Υ andΥ) and Substituting the (27) and (32) into (21) shows that as r → ∞.
To completely determine the form of the W r we must also determine the value of γ 1 . This requires matching the late-order expression given in (21) to the leading-order behaviour in the neighbourhood of the singularity. This procedure is described in Appendix B, and shows that From the results obtained in Appendix B, we find that as r → ∞. Hence, (35) reveals that the asymptotic series (8) may be expressed as the sum of two asymptotic series in even and odd powers of . In the following analysis, we apply optimal truncation methods to the asymptotic series, (8) and show that the optimal truncation error is proportional to U 3 e −χ3/ .

Stokes smoothing
In order to determine the behaviour of the exponentially small contributions in the neighbourhood of the Stokes curve we optimally truncate (8). We truncate the asymptotic series at the least even term by writing where N opt is the optimal truncation point and R N is the optimally-truncated error. As the analysis is technical, we will summarize the key results in this section with the details provided in Appendix C.
The leading order behaviour of the remainder term in the small limit can be shown to take the form where S 3 is the Stokes switching parameter which varies rapidly between constant values in the neighbourhood of Stokes curves. For algebraic convenience, we let T (s) represent the finite sum in (36) and use the shift notation denoted by F = F (s), F = F (s + ), F = F (s − ). By substituting the truncated series, (36), into the governing equation (5) we obtain where the terms neglected are of order O( 2N +1 W 2N +1 ) and quadratic in R N . Terms of these sizes are negligible compared to the terms kept in (38) as → 0.
Recall that Stokes switching occurs across Stokes curves, which are characterized by Im(χ) = 0 and Re(χ) > 0. We define the value of the Stokes multiplier, S 3 , to be In Appendix C we show that the Stokes multiplier, S 3 , changes in value by as Stokes curves are crossed and H is the function defined by H(χ 3 ) = (1 − 4W 3 0,3 ) 1/2 /χ 3 . Let S j denote the Stokes multiplier associated with χ j to be for j = 1, 2, 3 and C j is an arbitrary constant. Then the asymptotic power series expansion of Type A solutions of (5) with leading order behaviour, W 0,3 , up to exponentially small corrections is given by as → 0. By using the symmetry given by (7) we find that the other Type A solutions are given by as → 0, where S j is given by (40), λ 1 = e −2iπ/3 , λ 2 = e 2iπ/3 and In particular, the values of Λ j are where Λ 3 is given by (66). Furthermore, the complete asymptotic series expansions of Type A solutions of (5) have Stokes multipliers which contain a free parameter in the form of C j . This will be further explained in Section 3.3.
We have successfully determined a family of asymptotic solutions of (5) which contains exponentially small error. These exponentially small terms exhibit Stokes switching and therefore the expressions (41)-(43) describe the asymptotic behaviour in certain regions of the complex plane. The regions of validity for Type A solutions will be determined in Section 3.3.

Stokes structure
In this section we determine the Stokes structure of Type A solutions in both the complex s and xplanes. As demonstrated in Section 3.2, we found that the exponential contributions present in the series expansions (41), (42) and (43) are proportional to exp(−χ j / ). These terms are exponentially small when Re(χ j ) > 0 and exponentially large when Re(χ j ) < 0. Hence, by considering (27) we may determine the behaviour of these terms and the location of the (anti-) Stokes curves. Once the Stokes structure has been determined in the complex s-plane, we may reverse the scalings given by (4) in order to determine the Stokes structure in the complex x-plane.
We first illustrate and explain the Stokes structure and Stokes switching behaviour of (41) in the domain D 0 , described by (11). The upper and lower boundaries of D 0 are described by the curves Im(s) = ±π and are denoted by the dot-dashed curves in Figure 3.1a. In this figure we also see that there are two Stokes curves (red curves) and three anti-Stokes curves (dashed blue curves) emanating from the singularity s 0,3 . The Stokes curves extending towards the upper boundary of D 0 switches the exponential contribution associated with χ 3 as Re(χ 3 ) > 0. This Stokes curve is denoted by in Figure 3.1b. While the Stokes curve extending towards the lower boundary of D 0 does not switch any exponential contributions as Re(χ 3 ) < 0. Additionally, there is a branch cut (zig-zag curve) of χ 3 located along the negative real s axis emanating from s 0,3 . Using this knowledge, we can determine the switching behaviour as Stokes curves are crossed. Additionally, we observe in Figure 3.1a that the (anti-) Stokes curves and the branch cut separate D 0 into six regions.
We now determine regions in D 0 in which the asymptotic behaviour of (5) is described by the power series expansion (41), referred to as regions of validity. From Figure 3.1a we observe that the exponential contribution associated with χ 3 is exponentially small in the neighbourhood of the upper Stokes curves since Re(χ 3 ) > 0, and therefore the presence of exp(−χ 3 / ) does not affect the dominance of the leading order behaviour in (41). Hence, the value of the Stokes multiplier, S 3 , in the neighbourhood of the upper Stokes curve may be freely specified, and will therefore contain a free parameter hidden beyond-all-orders. The values of S 3 in the regions of D 0 is illustrated in Figure 3.2a.
The remainder term associated with χ 3 will exhibit Stokes switching and therefore the value of S 3 varies as it crosses a Stokes curve, say, from state 1 to state 2. Using the naming convention Re(s)   Figure 3.1a illustrates the behaviour of χ3 as Stokes and anti-Stokes curves are crossed. Figure 3.1b illustrates regions of D0 in which the exponential contribution associated with χ3 is exponentially large or small. The legend of this figure will apply to all subsequent figures unless stated otherwise.   (40) and (39), as Stokes curves are crossed. Figure 3.2b illustrates the regions of validity of the asymptotic solution described by (41). The dark and light grey shaded regions denote regions in which the exponential contribution associated with χ3 is large and small respectively. Hence, the regions of validity are those which are shaded in light grey.     Figure 3.1b extended to the adjacent domains D1 and D−1 as described by (12). In fact, the Stokes structure in the k th -adjacent domains are identical to those in D0 since W0 is 2πi-periodic. Each adjacent domain contributes an exponential contribution to the asymptotic solution described by (41). From figure 3.4a, we see that the presence of the exponential contribution from D−1 does not affect the dominance of (41) in D0 and hence its associated Stokes multiplier may be freely specified. However, the presence of the exponential contribution from D1 dominates the asymptotic solution, (41). In order for (41) to remain valid in D0, we require the value of the Stokes multiplier associated with exponential contribution in D1 must be chosen to be equal to zero. Figure 3.4b illustrates the regions where each exponential contribution are present.

Re(s)
described by (39), we denote these states by S − 3 and S + 3 respectively. If we assume that the value of S 3 is nonzero on either side of the upper Stokes curve, then we conclude that the exponentially small contribution associated with χ 3 is present in the regions bounded by the upper anti-Stokes curve, the anti-Stokes curve emanating from the singularity along the positive real s-axis and the curve Im(s) = π. Furthermore, the exponential contribution associated with χ 3 is also exponentially small in the region bounded by the branch cut and the lower anti-Stokes curve. The regions of validity of the asymptotic solution described by (41) are illustrated in Figure 3.2b.
However, for special choices of the free parameter hidden beyond-all-orders, we can obtain asymptotic solutions with an extended range of validity in D 0 . If we specify the value of S − 3 to be equal to zero, then the exponential contribution associated with χ 3 is no longer present in regions where it is normally exponentially large. In this case, the region of validity is extended by two additional adjacent sectorial regions in D 0 as illustrated in Figure 3.3b. We note that the case where S + 3 = 0 is specified can also give Type A solutions with an extended region of validity. However, this only extends the regions of validity of (41) by one additional sectorial region. In both cases the value of S 3 is specified, and therefore the asymptotic solution described by (41) is uniquely determined; we call these special Type A asymptotic solutions. Figure 3.4 illustrates the Stokes structure in the adjacent domains D 1 and D −1 as described by (12). Due to the 2πi-periodic nature of W 0 , the Stokes structure is also 2πi-periodic as shown in Figure 3 Figure 3.1, which is due to the symmetry described by (7). Figures 3.5a and 3.5b illustrate the regions of validity for the general and special asymptotic solution described by (42) . Figures 3.5c and 3.5d illustrate the regions of validity for the general and special asymptotic solution described by (43). D k do affect the asymptotic behaviour in D 0 . In order for the asymptotic solution (41) to correctly describe the solution behaviour in D 0 , the value of the Stokes multipliers must be specified such that the are not present in D 0 . The presence of the exponential contributions in the domains D −1 , D 0 and D 1 is illustrated in Figure 3.4b.

Re(x)
Im(x) (a) Stokes structure of (41) in the x-plane.

Re(x)
Im(x) (b) Regions of validity of (41) in the x-plane.

Re(x)
Im(x) (c) Extended regions of validity of (41) in the x-plane.

Re(x)
Im(x) (d) Stokes structure of (42) in the x-plane.

Re(x)
Im(x) (e) Regions of validity of (42) in the x-plane.

Re(x)
Im(x) (f) Extended regions of validity of (42) in the x-plane.

Re(x)
Im(x) (g) Stokes structure of (43) in the x-plane.

Re(x)
Im(x) (h) Regions of validity of (43) in the x-plane.

Re(x)
Im(x) (i) Extended regions of validity of (43) in the x-plane.  (44), the Stokes and anti-Stokes curves are now described by q-spirals in the x-plane. The branch cuts of the singulants, χ, in the x-plane extend from the singularities and terminate at the origin. As we are using the leading order term of the inverse transformation in the limit → 0, the boundaries of D0 are mapped to the negative real x axis, which corresponds to a logarithmic branch cut. Figure 3.6b illustrates the region of validity (light grey shaded regions) of (41) which contain a free parameter hidden beyond-all-orders. The regions of validity of (41) for which S − 3 is uniquely specified are illustrated in Figure 3.6c. The corresponding Stokes structure and regions of validity for the remaining two Type A solutions are illustrated in Figures 3.6d-3.6i.
The corresponding analysis of the Stokes structure and switching behaviour of the exponential contributions associated with the asymptotic solutions (42) and (43) may be obtained by using the symmetry (7). Consequently, the Stokes structure associated with χ 1 and χ 2 are vertical translations of the χ 3 -Stokes structure by ∓2iπ, respectively. The corresponding figures are contained in Figure 3.5.
As the Stokes structure and switching behaviour of the exponential contributions have been determined in the domain D 0 , we may finally determine the Stokes structure in the original complex x-plane. In order to determine the Stokes structure in the x-plane we reverse the scaling transformations. In particular, from the scalings given by (6) we find that the leading order mapping from s to x is given by as → 0. We illustrate the Stokes structures of Type A solutions for the choice of q = 1 + 0.2i. Using (44), the singulants, χ j (s), can be written as a function of x. We then compute the Stokes structure in the complex x-plane using Matlab. The corresponding Stokes structure for the asymptotic solution described by (41) in the complex x-plane is illustrated in Figure 3.6a. In these figures, the (anti-) Stokes curves and branch cuts in the complex x-plane follow the convention illustrated in Figure 3.1.
In particular, the complex x-plane contains a logarithmic branch cut (dot-dashed curve) along the negative real x axis as a result of using the leading order term of the inverse transformation described by (44). In fact, the boundaries of D 0 are mapped to this branch cut in the complex x-plane. The inverse transformation maps the Stokes and anti-Stokes curves in the s-plane to qspirals in the complex x-plane. From Figure 3.6, we find the Stokes and anti-Stokes curve separate the complex x-plane into sectorial regions bounded by arcs of spirals.

Type B Asymptotics
In this section we investigate Type B solutions of (5). These are the asymptotic solutions of (5), which are described by W 0,4 to leading order as → 0. The analysis involved in the subsequent sections is nearly identical to Sections 2 and 3. Hence, we will omit the details and only provide the key results.
Type B solutions may be expanded as a power series in of the form as → 0. Following the analysis for Type A solutions, it can be shown that the behaviour of y r is also described by a factorial-over-power form. The behaviour of y r is therefore described by as r → ∞. Recall that the main difference between Type A and Type B solutions is that Type B solutions are singular at three distinct points in D 0 rather than one. This feature will therefore be encoded in the calculation of the singulant function, η.
Since the leading order behaviour W 0,4 is singular at s 0,1 , s 0,2 and s 0,3 , we obtain three singulant contributions, which we denote by η j (s). This feature can be deduced from equation (17), which shows that W 0,4 is expressible as the sum of W 0,1 , W 0,2 and W 0,3 . Following the analysis in Section 3.1, η j (s) is given by for j = 1, 2, 3 and where σ is given in (26) with W 0,3 replaced by W 0,4 . Hence we obtain three contributions for η. Using the results for the late-order terms found in Section 3.1, we find that the late-order terms of (45) is given by as r → ∞, and where Y j (s) are the prefactor terms associated with η j (s). Hence the asymptotic expansion of Type B solutions of (5) is given by as → 0. As there are three distinct singulant terms in (48) there will be three subdominant exponentials present (after optimal truncation) and hence Type B solutions will also display Stokes behaviour. By applying the Stokes smoothing technique demonstrated in Section 3.2 to (48), the expression which captures the Stokes behaviour of the subdominant exponential correction terms is given by as → 0, and where N opt is the optimal truncation point. In particular, the Type B prefactor terms satisfy equation (32) with W 0 replaced by W 0,4 . Similarly, the Stokes multipliersŜ(s) satisfy (76) with W 0 and χ replaced by W 0,4 and η respectively. In view of the formula (17), the leading order behaviour of the Type B solution is a composition of the leading order behaviours of Type A solutions. Consequently, the Stokes behaviour present in this solution will be more complicated as the Stokes curves emanate from more than one point as this allows the possibility of interaction effects. In order to determine the Stokes structure of Type B solutions, we analyze the singulant (46).

Stokes Structure
The Stokes structure of Type B asymptotic solutions in D 0 is illustrated in Figure 4.1. In Figure  4.1 we see that there are three Stokes and two anti-Stokes curves curves emanating from each of the singularities, s 0,j , for j = 1, 2, 3. The Stokes structure for Type B solutions is more complicated as there are Stokes curves which cross into the branch cuts (zig-zag curves) of η j as illustrated in Figure 4.1. As these Stokes curves continue onto another Riemann sheet of η j , they may be subject to possible interaction effects from singularities originating in these Riemann sheets. However, we observe from   (5) in the domain D0. Regions I-IV denote regions in which the exponential contributions associated with ηj, which is given by (46), are exponentially small and therefore denote the regions of validity for (49).
in Figure 4.1. We first note that the Stokes curves emanating from the singularities s 0,1 , s 0,2 and s 0,3 asymptote to infinity as Re(s) → ∞ as illustrated in Figure 4.1a. In Figure 4.1a, we find that region I is the region bounded by the upper boundary of D 0 , the upper anti-Stokes curve emanating from s 0,1 and the Stokes curve emanating from s 0,2 , which is labelled by ; region II is the region bounded by the upper anti-Stokes curve emanating from s 0,1 , and the Stokes curves labelled by and ; region III is the region bounded by the lower anti-Stokes curve emanating from s 0,2 , and the Stokes curves labelled by and ; and finally region IV is the region bounded by the lower boundary of D 0 , the lower anti-Stokes curve emanating from s 0,2 and the Stokes curve emanating from s 0,1 , which is labelled by .
Furthermore, Re(η j ) is positive in each of these four regions and hence the exponential contributions associated with η j are exponentially small there. This is illustrated by the light grey shaded regions in Figure 4.1b. We therefore restrict our analysis to regions I-IV as these are the regions in which the dominant asymptotic behaviour is described by (49). Consequently, the regions of validity of Type B solutions are the regions bounded by the upper and lower anti-Stokes curves emanating from the singularities, s 0,1 and s 0,2 , respectively, and the boundaries of D 0 containing the positive real s axis. This is the union of regions I-IV and is illustrated in Figure 4.1b.
Hence, the Stokes curve separating regions I and II switches on the exponential contribution associated with η 2 , the Stokes curve separating regions II and III switches on the exponential contribution associated with η 3 and the Stokes curve separating regions III and IV switches on the exponential contribution associated with η 1 . To denote the Stokes switching behaviour of these subdominant exponentials, the Stokes curves are labelled by , and .
In regions I-IV, the presence of the exponential contributions associated with η j are exponentially small since Re(η j ) > 0, and therefore do not affect the dominance of the leading order behaviour in (49). Hence, the values ofŜ j may be freely specified in these regions and therefore the asymptotic solution described by (49) contains free parameters hidden beyond-all-orders. In Section 3.3 we were able to obtain special asymptotic solutions by uniquely specifying the free parameters present in the asymptotic expansion. However, this cannot be done for Type B solutions in the same way because of the possible interaction effects of singularities originating from the different Riemann sheets of η j .
Re  Under the inverse transformation (44) the Stokes and anti-Stokes curves are described by q-spirals in the x-plane. The exponential contributions associated with η1, η2 and η3 are switched across the Stokes curves labelled by , and , respectively as illustrated in Figure 4.2b. In Figure 4.2c the regions shaded in blue denote regions in which the asymptotic solution described by (49) is valid. These are the regions in which the exponential contributions present in (49) are exponentially small. The regions shaded in red are regions in which these exponential contributions are exponentially large and hence regions in which the asymptotic behaviour is not described by (49).
It may be possible to find special Type B solutions by considering the behaviour of the expo-nential contributions on the different Riemann sheets of η j and how they interact with those on the principal Riemann sheet. As this is beyond the scope of this study, we restrict our analysis to regions I-IV and therefore only obtain asymptotic solutions described by (49) which contain free parameters hidden beyond-all-orders. Following the analysis in Section 3.3 we obtain the Stokes structure in the original x-plane by applying the inverse transformation given by (44). As in Section 3.3 we demonstrate the Stokes structure for the value of q = 1 + 0.2i. Using Matlab, the Stokes structure of Type B solutions is illustrated in Figure 4.2. Figure 4.2 shows the corresponding regions I-IV in the complex x-plane. In this figure, the Stokes and anti-Stokes curves are denoted by the solid red and dashed blue curves respectively. The branch cuts of η j are depicted as the zig-zag curves, which connect the singularities to the origin in the complex x-plane. Furthermore, the dot-dashed curve denotes the logarithmic branch cut defined by the reverse transformation (44) for the choice of q = 1 + 0.2i.
Recall that the inverse transformation maps the Stokes and anti-Stokes curves in the complex s-plane to q-spirals in the complex x-plane. In particular, we see in Figure 4.2 that the Stokes curves of Type B solutions extend to infinity in the complex x-plane. This is due to the fact that the Stokes curves in the complex s-plane extend to infinity as Re(s) → ∞ as shown in Figure 4.1a. This was not case for the Stokes curves of Type A solutions in Section 3.3. Instead the Stokes curves of Type A solutions emanate from the singularities and approach the logarithmic branch cuts and enter a different Riemann sheet of the inverse transformation; this is illustrated in Figures  3.6a, 3.6d and 3.6g.
We have therefore determine the regions of validity for Type B solutions of q-P I , (5), in the complex x-plane. Type B solutions are described by the asymptotic power series expansion (49) as → 0, and contain free parameters hidden beyond-all-orders. Furthermore, we have also calculated the Stokes behaviour present within these asymptotic solution, (49), which allowed us to determine regions in which this asymptotic description is valid.

Connection between Type A and Type B solutions and the nonzero and vanishing asymptotic solutions of q-Painlevé I
In this section we establish a connection between both Type A and B solutions found in this study to the nonzero and vanishing asymptotic solutions of q-Painlevé I respectively. We recall that the nonzero and vanishing asymptotic solutions were first found by Joshi in [28]. In particular, equation (1) admits solutions with the following behaviour as |x| → ∞, where w nv and w v are the nonzero and vanishing asymptotic solutions and ω 3 = 1.
Recall that the four possible solutions for W 0 (s) are denoted by W 0,j (s), which are defined by (15)- (16). In our investigation, we applied the scalings given in (4), in which the variable x has the behaviour described by (6) as → 0. The analysis for both Type A and B solutions are valid in the limit → 0, which was shown to be equivalent to the double limit |q| → 1 and n → ∞. Under the additional limit, s → +∞, we see that the behaviour of x in (6) approaches infinity. Therefore, the limits → 0 and s → +∞ are equivalent to the limit |x| → ∞.
We now study the behaviour of W 0,j (s) under the additional limit s → +∞. By applying the limit s → +∞ to the terms appearing in (13) and (14) we find from equations (15) and (16) that The limiting behaviour of Type A solutions under the limit s → +∞ are therefore described by cube roots of unity. That is W 0,j ∼ ω 3 , as → 0 and s → +∞ for j = 1, 2, 3. However, we also find from (51) that Type B solutions vanish in the limit s → +∞. We therefore find that Type A solutions of (5) tend to the nonzero asymptotic behaviour solutions, w nv , found by Joshi [28] under the additional limit s → +∞.
In order to calculate the leading order behaviour of W 0,4 (s) as s → +∞, we need to keep terms up to order O(e −3s ). Be carefully tracking such terms we find that the behaviour of Type B solutions are given by as s → +∞. From (6) we find that (52) is equivalent to the behaviour 1/x as |x| → ∞ and hence Type B solutions correspond to the vanishing asymptotic behaviour, w v , of (1).

Numerical computation for q-Painlevé I
In this section we give a numerical example of (3) with the parameter choice of q = 1 + 0.2i and where x = x 0 q n with x 0 = 1. Given two initial conditions, w 0 and w 1 , a sequence of solutions of (3) may be obtained by repeated iteration. In general, only a certain choice of initial conditions will give a solution of (3) which tends to the asymptotic behaviour of interest. We follow the numerical method demonstrated by [31], originally based on the works of [29] to find appropriate initial conditions which tend to Type A solutions of (3).    Figure 5.1a and 5.1b show that real and imaginary part of w n converges to −1/2 and √ 3/2 respectively for large n, which is precisely described by the leading order term of w nv in (50).

Conclusions
In this study, we extended the exponential asymptotic methods used in [31,32] to compute and investigate Stokes behaviour present in the asymptotic solutions of q-P I in the double limit |q| → 1 and n → ∞. In order to investigate the solution behaviour of q-difference equations we rescaled the variables in the problem such that the double limit |q| → 1 and n → ∞ was equivalent to the limit → 0. We found two types of solutions for q-P I , which we call Type A and Type B solutions. Type A solutions of q-P I are described asymptotically by either (41), (42) or (43), and Type B solutions are described by (49). The asymptotic descriptions obtained in this analysis are given as a sum of a truncated asymptotic power series and an exponentially subdominant correction term. We then determined the Stokes structure and used this information to deduce the regions of the complex plane in which these asymptotic solutions are valid.
In Section 3 we first considered the asymptotic solutions of (5) described by Type A solutions. Using exponential asymptotic methods, we determined the form of the subdominant exponential contribution present in the asymptotic solutions, which were found to be defined by one free Stokes switching parameter. From this behaviour, we deduced the associated Stokes structure, illustrated in Figures 3.1a and 3.5. By considering the Stokes switching behaviour of these subdominant exponentials, we found that the dominant asymptotic behaviour is described by either (41), (42) or (43) in a region in the complex s-plane containing the positive real axis. Furthermore, we found that it is possible to select the Stokes parameters so that the exponential contribution is absent in the regions where it would normally be large. Consequently, the regions of validity for the special Type A solutions are larger than the regions of validity for generic Type A solutions as illustrated in Figures 3.3b, 3.5b and 3.5d.
In Section 4, we considered the equivalent analysis for Type B solutions of (5). Compared to Type A solutions, Type B solutions are those which are singular at three points rather than one. Although exponential asymptotic methods may be used again to determine the form of the exponential small contributions present in Type B solutions, we noted that the analysis is near identical except for the determination of the singulant of this problem, η. The main difference is due to the fact that Type B solutions are singular at three points rather than one, and the remaining analysis was therefore identical as for Type A solutions. Type B asymptotic solutions are given as a sum of a truncated asymptotic power series and three exponentially subdominant correction terms as described in (49). The Stokes structure for Type B solutions, illustrated in Figure 4.1, is significantly more complicated than the Stokes structure of Type A solutions. In order to describe the Stokes switching behaviour in the domain D 0 we must understand how the asymptotic solution (49) interacts with solutions on different Riemann sheets. However, we restrict ourselves to regions I-IV as these regions are free of possible interaction effects with singularities originating from the different Riemann sheets. Furthermore, all three exponential contributions present in Type B solutions are exponentially subdominant in regions I-IV and therefore represent the regions of validity of (49). Consequently, the asymptotic solutions described by (49) contain one free parameter defined by the Stokes multiplier.
By reversing the rescaling transformations we were then able to obtain the corresponding Stokes structure in the complex x-plane and found that the Stokes and anti-Stokes curves are described by q-spirals in the complex x-plane. As a result, the regions of validity are no longer described by traditional sectors bounded by rays but sectorial regions bounded by arcs of spirals. Consequently, this methodology is applicable to other q-Painlevé equations or more generally nonlinear q-difference equations in order to obtain asymptotic solutions which display Stokes phenomena.
Finally, Section 5 demonstrated that Type A and B solutions are related to the nonzero asymptotic and vanishing asymptotic solutions found by Joshi [28]. If the additional limit s → +∞ is also taken, then we find that Type A solutions correspond to the nonzero asymptotic solutions of q-P I while Type B solutions correspond to the quicksilver solutions of q-P I [28].
Appendix A Calculating the late-order terms near the singularity In order to determine the value of γ 1 in (21) we calculate the local behaviour of χ 3 and U 3 . Using (20) in equations (26) and (32) we can show that as s → s 0,3 . Furthermore, by calculating the expression for W 1 (using equation (9)) it can be shown that W 1 (s) = O (s − s 0,3 ) −1/2 and from (33) that as s → s 0,3 for some constant A 1 . By substituting (53) into the expression for W r as given by (34) we find that the local behaviour of of the late-order terms near the singularity is given by as s → s 0,3 . From (9), the dominant behaviour of W r is due to the term W r−2 /(4W 3 0 − 1). Hence, if W r−2 has a singularity of strength ν, then W r has one of strength ν + 5/2. In particular, since we know that W 0 is singular at s 0,3 of strength −1/2 then W 2r will also be singular at s 0,3 but with strength 5r/2 − 1/2. Moreover, since W 1 is also singular at s 0,3 of strength 1/2, then W 2r+1 has one of strength 5r/2 + 1/2.
with E 0 = 1. We then substitute (62) into (61) and match terms of O(ζ) in the limit |ζ| → ∞. Doing this, we obtain the following nonlinear recurrence relation for r ≥ 1. We can express (62) in terms of the outer variables by reversing the scalings given in (58). Doing this, we find that the inner expansion of the outer solution is given by as → 0. Recall that the outer expansion is given by the expression (8) as → 0, and where the behaviour of coefficients are given by (35). By matching the expansions (8) and (64) it follows that We then compute the first 1000 E r terms using the recurrence relation (63). Then, by using the formula (65) we find numerically that the approximate value of Λ 3 is to four significant figures. The approximate value for Λ 3 is shown in Figure B approximation for Λ appearing in the even late-order terms (55). As r increases, the approximation for Λ tends to the limiting value of −0.04364, which is denoted by the black dashed line.

Appendix C Stokes smoothing
To apply the exponential asymptotic method, we optimally-truncate the asymptotic series (8).
One particular way to calculate the optimal truncation point is to consider where the terms in the asymptotic series is at its smallest [8]. This heuristic is equivalent to the finding N such that 2N +2 W 2N +2 2N W 2N ∼ 1, in the limit N → ∞. By using the late-order form ansatz described by (21) we find that N ∼ |χ 3 |/(2 ) as → 0 (which is equivalent to the limit N → ∞). As this quantity may not necessarily be integer valued, we therefore choose κ ∈ [0, 1) such that is integer valued. By substituting the truncated series, (36), into the governing equation (5) we obtain T T 2 T + T 2 T R + 2T T T R + T T 2 R + · · · = T + R − 1 where the terms neglected are quadratic in R N . We can use the recurrence relation (9) to cancel terms of size o( 2N −1 ) in (68). In particular, equation (68) can be rewritten, after rearrangement, as as → 0 and where the terms neglected are of order O( 2N +1 W 2N +1 ) and terms quadratic in R N . Terms of these size are negligible compared to the terms kept in (69) as → 0. Away from the Stokes curve the inhomogeneous terms of equation (69) are negligible as → 0. We may therefore apply a WKB analysis to the homogeneous version of (69) by setting R N,hom (s) = α(s)e β(s)/ as → 0. Using the WKB method, we find that the leading order equation gives the solution β = −χ 3 . Continuing to the next order involves matching terms of O( R N ). By collecting terms of this size in the homogeneous version of (69) we obtain the equation − 2W 3 0,3 sinh(χ 3 ) α α − W 3 0,3 χ 3 cosh(χ 3 ) + 2W 2 0,3 W 0,3 sinh(χ 3 ) + 3 where we have substituted the fact that β = −χ 3 . Comparing equations (70) and (28) show that α satisfies to same differential equation as the prefactor U 3 and hence α ∝ U 3 . Hence, away from the Stokes curves the solution of (69) is given by R N,hom (s) ∼ U 3 (s)e −χ3(s)/ , as → 0. In order to determine the Stokes switching behaviour associated with the optimally-truncated error, we set R(s) = S 3 (s)R N,hom (s), where S 3 (s) is the Stokes multiplier. We substitute (71) into (69) and use equations (25) and (28) to cancel terms. Doing this we find that the Stokes multiplier satisfies as → 0. Noting the form of N , we introduce polar coordinates by setting χ 3 = ρe iθ where the fast and slow variables are θ and ρ, respectively. This transformation tells us that and (67) becomes N opt = (ρ/ + κ)/2. Under this change of variables equation (72) becomes as → 0, where we have used Stirling's approximation [1] of the Gamma function. For simplicity, we will let H(s(θ); ρ) = 1 − 4W 3 0,3 /χ 3 . From (74) we find that the right hand side is exponentially small everywhere except in the neighbourhood of θ = 0, which is exactly where the Stokes curve lies (where χ 3 is purely real and positive). We now rescale about the neighbourhood of the Stokes curve in order to study the switching behaviour of S 3 by setting θ = √ θ . Note that under this scaling, H(s(θ); ρ) ∼ H(|χ 3 |) as → 0, which is therefore independent of θ to leading order. Applying the scaling θ = √ θ to (74) gives as → 0. Integrating (75) we find that where C 3 is an arbitrary constant. Thus, as Stokes curves are crossed, the Stokes multiplier changes in value by ∆S 3 ∼2iπ √ H(|χ 3 |), and therefore ∆R N ∼ 2iπ √ H(|χ 3 |)U 3 (s)e −χ3(s)/ , as → 0.