On the origins of Riemann–Hilbert problems in mathematics

This article is firstly a historic review of the theory of Riemann–Hilbert problems with particular emphasis placed on their original appearance in the context of Hilbert’s 21st problem and Plemelj’s work associated with it. The secondary purpose of this note is to invite a new generation of mathematicians to the fascinating world of Riemann–Hilbert techniques and their modern appearances in nonlinear mathematical physics. We set out to achieve this goal with six examples, including a new proof of the integro-differential Painlevé-II formula of Amir et al (2011 Commun. Pure Appl. Math. 64 466–537) that enters in the description of the Kardar–Parisi–Zhang crossover distribution. Parts of this text are based on the author’s Szegő prize lecture at the 15th International Symposium on Orthogonal Polynomials, Special Functions and Applications (OPSFA) in Hagenberg, Austria.

(Some figures may appear in colour only in the online journal)

Introduction
Riemann-Hilbert problems, in short RHPs, are useful in the analysis of orthogonal polynomials (OPs), special functions (SFs) and thus several applications (A) in mathematics or physics. Still, these modern appearances of RHPs are quite detached from the original meaning of the term RHP and it is our primary objective to pin down the occurrence of the first problem in mathematics that was coined an RHP. Once done, we will follow Plemelj's traces whose work on the original RHP gave rise to an analytic apparatus, the Riemann-Hilbert techniques or Riemann-Hilbert methods, which found various applications in the OPSFA realm. We will showcase and advertise the latter modern aspects of the theory with examples from mathematical physics, in particular examples from nonlinear wave and integrable systems theory, statistical mechanics, random matrix theory and integrable probability.

The original Riemann-Hilbert problem
In order to answer the question 'What is . . . a Riemann-Hilbert problem' one should go back to 1900, the year of the second international congress of mathematics. Held at the turn of the century in Paris, the ICM 1900 proved to be highly influential for the development of 20th century mathematics, in particular since Hilbert delivered his famous address in the section Bibliographie et Histoire. Enseignement et methodes. During his lecture he presented (because of time constraints) ten mathematical problems 'from the discussion of which an advancement of science may be expected', see [70, p 445]. The full, twenty-three long, list of Hilbert's problems appeared in German print even before the French proceedings of the congress. Scrolling through this list we eventually encounter problem #21, the main focus point of this section, see [69, pp 289-290]:
In this problem Hilbert asks the reader to show that there always exists a linear differential equation of the Fuchsian class, with given singular points and monodromic group. Before reviewing some of the necessary background material in section 3 below it will be important to point out the following peculiar wording in Hilbert's original formulation: the problem requires the production of n functions of the variable z, regular throughout the complex z plane except at the given singular points; at these points the functions may become infinite of only finite order. Thus, Hilbert does not use the word pole singularity-he refers to what is nowadays called a regular singularity, see [83, definition 16.2]-and he also never speaks of Fuchsian linear systems. These subtleties in the formulation of the problem played an important role as we will soon learn.
Right now we follow Anosov and Bolibrukh [5, p 7] and coin Hilbert's 21st problem the original RHP: Hilbert formulated it and he explicitly mentions Riemann in his wording as this problem is extremely close to Riemann's ideas 1 (namely the global construction of a function from given local analytic data, i.e. the singularity locations and associated monodromy). In short,

Riemann-Hilbert Problem 2.1 (The original RHP ([69], 1900)) . Proof of the existence of linear differential equations having a prescriped monodromic group.
This problem seems quite different from the modern RHPs discussed in section 6 below, yet it is this problem on which the modern theory is based upon. Thus, being so influential, it is only natural to ask if RHP 2.1 has been solved 120 years after Hilbert's address? Well, according to the Encyclopedia of Mathematics [51], the 'solution is negative or positive depending on how the problem is understood'. This answer points at a curious mathematical misunderstanding which persisted from 1908 till 1983 and which relates to the aforementioned atypicalities in Hilbert's wording. Indeed, we highlight two possible sources of confusion: • equation vs. system: The monodromy of a pth order scalar linear ODE with solution y matches the monodromy of the p × p linear ODE system which the vector (y, y , . . . , y p−1 ) solves. Still, at the time of the congress, see [123], it was known that Hilbert's question for scalar Fuchsian equations has a negative solution unless additional singularities are included. So by a pth order scalar Fuchsian equation, Hilbert likely meant a generic p × p Fuchsian system. • regular singularities vs Fuchsian singularities: Both notions exist for pth order scalar linear ODEs in the complex plane. Although up front different in appearance, these notions are in fact equivalent (this is a result by Fuchs from 1868, see [83, theorem 19.20]). However, this is no longer the case when one considers the same singularities for p × p linear ODE systems, see section 4.
In turn we could think of at least three different interpretations of RHP 2.1: given a monodromy group with encoded singularity locations, are we asked to realize it by (a) a Fuchisan linear pth order differential equation? (b) a linear p × p system having only regular singularities? (c) a Fuchsian system on the whole Riemann sphere CP 1 ?
The answer to (a) is negative as shown by Poincaré [123]: the number of parameters in a pth order linear Fuchsian equation with singularities a 1 , . . . , a n is less than the dimension of the space M of monodromy representations 2 . Interpretation (b) was positively solved by Plemelj in 1908 [121] and believed to settle (c) as well. This belief persisted for 75 years until Kohn [94] and Arnold and Il'yashenko [7] discovered a serious gap in Plemelj's argument. As it turned out, interpretation (c) is much more subtle and a negative answer to Hilbert's 21st problem in said context was eventually given by Bolibrukh in 1989, see [17,18]. We will now discuss in detail interpretations (b), (c) and highlight the pioneering work of Plemelj in this context. His contributions laid the foundation for several aspects of modern Riemann-Hilbert theory which we showcase through a series of OPSFA-type examples in sections 6, 8 and 9-however only section 9 contains original content. For an in-depth discussion of interpretation (a) we refer the interested reader to [5, chapter 7] or more recently [66, section 3].

Background terminology
We are considering p × p linear ODE systems of the form where A(z) ∈ C p×p is analytic in a disk D r (z 0 ) ⊂ C of radius r > 0 centred at some z 0 ∈ C. Then, by the Picard-Lindelöf theorem [83, theorem 1.1], any fundamental solution of (3.1) is analytic in D r (z 0 ). Now suppose A(z) is analytic on the punctured Riemann sphere S := CP 1 \{a 1 , . . . , a n } where a 1 , . . . , a n are n distinct points on CP 1 . Through the use of overlapping disks and the previous local existence theorem, any fundamental solution of (3.1) can be continued along any path in S. This statement is part of the monodromy theorem [10, theorem 67] which furthermore asserts that the continuation only depends on the homotopy class of the path. With this in mind, we now define the monodromy of (3.1) as follows: fix a 0 ∈ S distinct from a 1 , . . . , a n and continue any fundamental solution Ψ(z) of (3.1) along some γ ∈ π 1 (S, a 0 ) (the fundamental group of the surface S with base point a 0 ). The continuation yields Ψ(z), another fundamental solution of (3.1), which links to Ψ(z) via Ψ(z) = Ψ(z)G γ , for some invertible p × p matrix G γ , i.e. G γ ∈ GL(p, C). Strictly speaking, G γ only depends on the homotopy class of the loop γ. In turn, the map [γ] → G γ defines a representation We are left with two obvious ambiguities in our definitions: first, if we were to start from a different fundamental solution Ψ * (z) instead of Ψ(z), then Ψ * (z) = Ψ(z)H for some H ∈ GL(p, C) and the associated monodromy matrices change accordingly, Second, the dependence of π 1 (S, a 0 ) and thus G i = χ([γ i ]) on the base point a 0 : given that S is path connected we know that π 1 (S, a 0 ) and π 1 (S, a * 0 ) are isomorphic for any two base points a 0 , a * 0 different from a 1 , . . . , a n , thus under the representation (3.2) the associated monodromy matrices change also in the style of (3.3). Hence, summarizing the two ambiguities, the monodromy of (3.1) is defined up to conjugation equivalence and so an element of the space M := Hom (π 1 (S, a 0 ), GL(p, C))/ GL(p, C), of conjugacy classes of representations of π 1 (S, a 0 ). Alternatively, in the language of monodromy matrices, we have M ∼ = {(G 1 , . . . , G n ) : G 1 · . . . · G n = I}/ GL(p, C). (3.4) One notion remains: we say system (3.1) is Fuchsian if all its singular points {a 1 , . . . , a n } are first order poles of A(z), i.e. without loss of generality (thanks to a fractional linear map) we have Using the above terminology we are now able to rephrase RHP 2.1 in interpretation (c) as follows: of Fuchsian systems with fixed singularities a 1 , . . . , a n into M surjective?

Riemann-Hilbert Problem 3.1 (The original RHP, case (3)). Is the monodromy map
This problem concludes our content on the necessary background terminology of Hilbert's 21st problem. We will now review some of the mathematical works devoted to its solution and which were published between 1908 and 1989. These works introduced several key ideas and techniques that are nowadays used in the Riemann-Hilbert analysis of OPSFA-type problems.

Plemelj's contributions
In 1908, Plemelj [121] (upgraded to book form 55 years later in [122]) published a solution of RHP 3.1 that was widely accepted until the early 1980s. In his work, he reduced RHP 3.1 to a Hilbert boundary value problem 3 in the theory of singular integral equations and the idea goes as follows: join all singularities a 1 , . . . , a n ∈ C by a simple closed oriented contour Γ as indicated in figure 2 below. which, thanks to the cyclic constraint (3.4), satisfies G(z) = I for z ∈ [a n , a 1 ). With Ω + denoting the region in C bounded by Γ and Ω − the complement of the closure Ω + in CP 1 , we then consider the following boundary value problem:

Hilbert Boundary Value Problem 4.1 ([121, pp 213, 214, 229])
. Find all vector-valued functions X = X(z) ∈ C 1×p such that (a) X(z) is analytic in Ω ± \{∞} and extends continuously from either side to the punctured contour Γ\{a 1 , . . . , a n }. satisfy the boundary condition X + (z) = X − (z)G(z), compare figure 3 below. 3 In standard textbooks on singular integral equations, see [112, section 39], a Riemann-Hilbert problem, named after the original works [71,72,125], generally refers to the problem of constructing a function which is analytic in a domain Ω ⊂ C, continuous on the closure Ω and with prescribed boundary values on ∂Ω. We will not follow this tradition but instead use the term Hilbert boundary value problem like in [138, section 4] as to not confuse ourselves with the original RHPs 2.1 and 3.1. (c) X(z) is of finite degree at z = ∞, that is With problem 4.1 as starting point, Plemelj first proved its solvability through an application of Fredholm's theory of integral equations, see [61,119], at the time a novel analytic tool. We will lay out some of his steps below subject to the following two temporary assumptions: X − (z), z ∈ Γ is the boundary value of some X(z) ∈ C 1×p , analytic in Ω − \{∞}, Hölder continuous on (Ω − \{∞}) ∪ Γ and of finite degree at z = ∞, in the sense of (4.1). Next, Plemelj noticed that both constraints in (4.2) are equivalent (again by Cauchy's theorem and the Plemelj-Sokhotski formulae) to the principal value integral equations, and from these one obtains in turn the quasi-regular integral equation with the Hölder continuous matrix-valued kernel function

Remark 4.4.
In more abstract and general terms, the left-hand side in (4.4) defines a singular Fredholm operator of index zero which acts on Hölder continuous functions defined on Γ, see [112, section 45], and we assume that K(z, λ) satisfies a Hölder condition on Γ in both variables. Plemelj did not work with such singular Fredholm equations in [121] since he managed to transform the piecewise constant jumps in problem 4.1 to differentiable ones. In turn, his [121, (4)] is an ordinary Fredholm equation of the second kind which can be studied by Fredholm's 1903 theory. We will discuss a more generally applicable reduction of piecewise Hölder continuous jumps to Hölder continuous ones and en route lift assumption 4.3 on G(z). In fact, the general theory of singular Fredholm integral equations of the type (4.4) was systematically developed some 30 years after Plemelj's paper with prominent contributions by Giraud [64], Muskhelishvili-Vekua [113] and Gakhov [62].
We now make (4.4) the starting point for the solvability analysis of problem 4.1. Precisely, given that any continuous solution X − (z), z ∈ Γ of (4.4) will be automatically Hölder continuous as a consequence of assumption 4.3, we answer the following two questions: Q1 Under which conditions is (4.4) solvable in the space of continuous functions on Γ? Q2 Does each continuous solution of (4.4) yield a solution of problem 4.1?
Regarding Q2, we recall that every continuous solution X − (z), z ∈ Γ of (4.4) will produce a solution of problem 4.1 with assumption 4.3 in place if and only if X − (z), z ∈ Γ solves both equations in (4.3). Now rephrase (4.3) in terms of the function which is analytic on Ω ± , Hölder continuous on Ω ± ∪ Γ and vanishes at z = ∞. In fact, Furthermore, our central integral equation (4.4) is equivalent to the jump condition which motivates the introduction of the below, after Plemelj the accompanying, Hilbert boundary value problem: The special homogeneous problem 4.5 is useful since it allows us to summarize our previous chain (4.6) as the following exclusive alternative: let X − (λ), λ ∈ Γ be a continuous solution of (4.4), • If X(z) defined in terms of said X − (λ), λ ∈ Γ in (4.5) is identically zero, then X − (λ) solves (4.3) and thus produces a solution of the initial problem 4.1.
non-trivial solution of the accompanying problem 4.5.
Most importantly we can now formulate Plemelj's first criterion as response to Q2 above: Evidently, lemma 4.6 does not guarantee solvability of (4.4). Indeed, in order to answer Q1, Plemelj used Fredholm's theory of integral equations 5 . First, consider two additional homogeneous Hilbert boundary value problems: the associated Hilbert boundary value problem which consists in finding Y = Y(z) ∈ C 1×p that is analytic in Ω ± , extends Hölder continuously up to Γ, vanishes at infinity and instead of the jump condition in problem 4.5 satisfies with the matrix transpose G (z) of G(z). Moreover, the problem which accompanies the associated problem (4.7) and which seeksŶ =Ŷ(z) ∈ C 1×p as in problem 4.5 but with jump constraintŶ (4.8) 5 Fredholm's theorems for singular integral equations of the form (4.4) were proven in [115] and are called Noether's theorems. If however the index of the underlying Fredholm operator is zero, compare remark 4.4, then those theorems are exactly the same as for standard Fredholm integral equations of the second kind, see [61,119].
Repeating the logic that took us from problem 4.1 (with assumptions 4.2 and 4.3) to equation (4.4) we easily see that (4.8) leads to the integral equation (4.9) and which is the adjoint equation of (4.4). We can now derive Plemelj's second criterion: Proof. Problem (4.8) accompanies the associated problem (4.7), or equivalently, (4.7) is the accompanying problem of (4.8) (as we invert the jump matrix and homogenize at infinity in the accompanying problems). But by assumption, (4.7) has no non-trivial solutions, so by lemma 4.6 every continuous solutionŶ + (z), z ∈ Γ of the adjoint equation (4.9) yields a solution of the Hilbert boundary value problem (4.8). LetŶ + (z), z ∈ Γ denote such a non-trivial continuous solution of (4.9), 6 soŶ + (z), z ∈ Γ is the boundary value of some analytic functionŶ(z), z ∈ Ω + . However, by [115, (7)], the inhomogeneous equation (4.4) is continuously solvable if and only if Γ γ(z)Ŷ + (z) dz = 0, and which is guaranteed by our last conclusion and the fact that γ ∈ C 1×p [z] (any entire γ(z) would do). This concludes the proof of the lemma.
To summarize, lemmas 4.6 and 4.7 guarantee solvability of our initial problem 4.1 (with assumption 4.3 in place, but no longer assumption 4.2), provided the accompanying problem 4.5 and the associated problem (4.7) are only trivially solvable. Although these two conditions seem peculiar they are in fact easily verifiable from our assumption that γ(z) is polynomial and thus our solutions of problem 4.1 meromorphic at z = ∞. Indeed, we first recall, see [115, section 2], that the homogeneous version of (4.4) (or any of the other homogeneous equations in this section) has finitely many, say s ∈ Z 1 , linearly independent continuous solutions 7 . Then Proof. If X(z) has a zero of order k ∈ Z 1 at z = ∞, then all solve the homogeneous version of problem 4.5 (with assumption 4.3) and thus (4.4) with γ ≡ 0. Since these solutions are linearly independent we must have k s by our previous discussion.
So we can always find an integer r ∈ Z 0 such that neither the accompanying problem 4.1 nor the associated problem (4.7) admit solutions with a zero of order greater than r at z = ∞. Hence, if X(z) solves problem 4.1 where γ ∈ C 1×p [z] satisfies deg(γ) r, then (a) The accompanying problem has jump X • so if this problem has a non-trivial solution (that vanishes at infinity) then solves problem 4.5 and has a zero of order greater than r at z = ∞, contradicting our initial choice for r.
so if this problem has a non-trivial solution (that vanishes at infinity) then solves the associated problem (4.7) and vanishes to order greater than r at z = ∞. Again a contradiction to our choice of r.
The X • -Hilbert boundary value problem is thus solvable by lemmas 4.6 and 4.7 and by reversing (4.10) we arrive at with some constants γ 1 , . . . , γ m ∈ C and where {X i (z)} m j=1 are linearly independent particular solutions of problem 4.1 with and {X i (z)} m i=p+1 (if there are any) have degree strictly less than r at z = ∞.
Proof. Clearly, if {X i } m j=1 are any particular solutions of the jump condition in problem 4.1, then any polynomial linear combination of them will also solve the same jump constraint. Thus, choosing r ∈ Z 1 sufficiently large as indicated in our previous discussion we can always ensure that problem 4.1 is solvable for the collection of constant polynomials γ i (z) = [δ i j ] p j=1 with i = 1, . . . , p. Reversing subsequently (4.10) and taking an arbitrary combination of the so obtained solutions yields the first half in (4.11) with the indicated asymptotic behaviour at z = ∞. The second sum in (4.11) stems from the complete system of linearly independent solutions of the homogeneous version of problem 4.1, modulo the inversion of (4.10).
Although theorem 4.9 solves problem 4.1, the implicit dependence of (4.11) on r ∈ Z 1 is a drawback, in particular since r affects X • (z) through (4.10) and it affects the parameter m. In order to bypass these deficiencies Plemelj managed then to construct a fundamental system of solutions to problem 4.1 such that every solution of the same problem is a polynomial linear combination of the form In order to present his argument, we first note that the system {X i } p i=1 in (4.11) is linearly independent over C[z] because of its asymptotic normalization. With r ∈ Z 1 as stated in theorem 4.9, we then observe that among all solutions of the form (4.11) some have the lowest degree at z = ∞ since a solution of problem 4.1 cannot vanish to all orders at z = ∞. Let F 1 (z) denote one of the solutions of lowest degree −ν 1 ∈ Z at z = ∞ and F 2 (z) one of the solutions of problem 4.1 of lowest degree −ν 2 −ν 1 at z = ∞ such that F 2 / ∈ span{F 1 }. 8 Next, we let F 3 (z) denote one of the solutions of problem 4.1 of lowest degree −ν 3 −ν 2 at z = ∞ such that F 3 / ∈ span{F 1 , F 2 } and note that this process may be continued until we reach . . , F k (z) with k < p have been constructed in this fashion then we can find yet another solution of problem 4.1 not contained in span{F 1 , . . . , F k }. For otherwise, the particular solutions in (4.11) would be contained in span{F 1 , . . . , F k } and thus be linearly dependent over C[z], contradicting our previous discussion. In summary, the above procedure leads to p fundamental solutions where {F i } p j=1 is a fundamental system of solutions, q i ∈ C[z] and the corresponding fundamental matrix (a) det Φ(z) = 0 for all z ∈ C including z ∈ Γ with the appropriate limiting values Φ ± (z).
Proof. The first statement of the theorem was already established in theorem 4.9 and (4.13) follows from (a) since by construction of the fundamental matrix, , z ∈ Γ and thus for any solution X(z) of problem 4.1 (with assumption 4.3), i.e. X(z)(Φ(z)) −1 is an entire function. Hence, by condition (c) in problem 4.1 we find that , so (4.13) follows and conversely it is clear that the right-hand side in (4.13) solves the Hilbert boundary value problem 4.1 for any q i ∈ C[z]. In order to establish (a) we first note that if X(z) is any solution of problem 4.1 of degree strictly less than −ν k ∈ Z at z = ∞ where k ∈ {2, . . . , p}, then necessarily (4.14) For if such X, i.e. of degree strictly less than −ν k at z = ∞, were not contained in span{F 1 , . . . , F k−1 }, then our previous selection process of the subsequent F k / ∈ span{F 1 , . . . , F k−1 } with lowest degree −ν k −ν k−1 would fail. Returning to (a) we now show that the expression with a i ∈ C not all zero, does not vanish at any point z ∈ C: if F(z 0 ) = 0 for some z 0 / ∈ Γ, then where X(z) solves problem 4.1. If a k denotes the last of the coefficients a 1 , . . . , a p which is non-zero, then the degree of X(z) at z = ∞ is strictly less than −ν k , so (4.14) must hold for the same X(z), i.e. we would have But this is a contradiction to the requirement F k / ∈ span{F 1 , . . . , F k−1 } and thus F(z) = 0 for all z ∈ C\Γ. If on the other hand F(z 0 ) = 0 for some z 0 ∈ Γ, i.e. either F + (z 0 ) = 0 or equivalently F − (z 0 ) = F + (z 0 )(G(z 0 )) −1 = 0, then define This C 1×p -valued function is analytic in Ω ± \{∞}, of finite degree at z = ∞ and extends Hölder continuously from either side to the punctured contour Γ\{z 0 }. Thus the above X(z) solves the singular integral equation (4.4) for all z ∈ Γ\{z 0 } and since F ± (z), z ∈ Γ are Hölder continuous with F ± (z 0 ) = 0 we conclude that X − (z), z ∈ Γ as well as X + (z) = X − (z)G(z), z ∈ Γ are Hölder continuous on all of Γ. At this point we simply repeat the previous logic for z 0 / ∈ Γ and conclude F(z) = 0 for all z ∈ Γ. All together, (4.15) does not vanish for any z ∈ C (choosing F ± (z) for z ∈ Γ) and thus det Φ(z) = 0 for all z ∈ C, as claimed. We are left with (b) and thus the point z = ∞ at which det Φ(z) might have a pole or may remain finite. However, if −ν i ∈ Z denotes the degree of F i (z) at z = ∞, then is certainly finite at z = ∞. If it were zero, then we could find a 1 , . . . , a p ∈ C not all zero such that With a k denoting the last of the coefficients a 1 , . . . , a p which is non-zero we then have i.e. the degree of X(z) at z = ∞ is strictly less than −ν k . But then necessarily (compare above) and in turn F k ∈ span{F 1 , . . . , F k−1 }, a contradiction. All together, This concludes the proof of the theorem.
Plemelj's central theorem 4.10 will be the key to solve RHP 2.1 in interpretation (b), provided we can dispose of assumption 4.3. In [121, pp 229-236] this was achieved by relying on the piecewise constant nature of G(z) in problem 4.1. Here, we shall outline a procedure of Vekua [138, chapter 2] which is more widely applicable and which employs multi-valued analytic functions defined in cut planes, equivalently discontinuous single-valued functions. First, we denote with G(a j ∓ 0), j = 1, . . . , n the left, resp. right limits at the point a j ∈ Γ as we approach a j along (a j−1 , a j ), resp. (a j , a j+1 ) on Γ, compare figure 2. These limits exist by our definition of G(z) and we call a σ ∈ {a 1 , . . . , a n } with G(a σ − 0) = G(a σ + 0) a point of discontinuity. Second, we use the eigenvalues λ σ 1 , . . . , λ σ p of the non-singular matrix G(a σ − 0)(G(a σ + 0)) −1 at a point of discontinuity a σ written as Third, we fix z 0 ∈ Ω + and let L σ denote the straight line connecting z 0 with a point of discontinuity a σ and extending to z = ∞. With this convention any branch of is single-valued and we have as well as We also need an arbitrary branch of the function (z − a σ ) γ σ j defined in C with a cut along We now come to the heart of Vekua's method: define the diagonal matrices for properly chosen M σ , N σ ∈ GL(p, C), will solve a modified problem 4.1 with jump matrix which is Hölder continuous on the whole segment (a σ−1 , a σ+1 ) including the point z = a σ .
Once a σ has been removed in this fashion all remaining points of discontinuity can be treated similarly and we eventually arrive at a Hilbert boundary value problem that satisfies assumption 4.3. So how do we choose M σ and N σ ? With (4.16) the jump condition on Γ\{a 1 , . . . , a n } in problem 4.1 gets replaced by and this motivates the requirements to ensure that the jump for T(z) is Hölder continuous on (a σ−1 , a σ+1 ). In fact, with (4.17) in place, the matrix (C σ (z)) −1 N σ G(z) − M σ will be Lipschitz continuous on (a σ−1 , a σ+1 ), see [112, section 6]. In order to make (4.17) more transparent, we use C σ (a σ − 0)(C σ (a σ + 0)) −1 = p j=1 λ σ j and rewrite (4.17) as If this is not the case, then one has to properly modify M σ , N σ in (4.16) and work with the Jordan normal form of G(a σ − 0)(G(a σ + 0)) −1 , see [138, pp 87-93] for details. Either way, a finite number of transformations of the type (4.16) lead to a Hilbert boundary value problem that is amenable to our previous analysis, i.e. a problem whose jump matrix is Hölder continuous along all of Γ 9 and which satisfies a normalization of the form (4.1). At this point we can come full circle and return to the original RHP 2.1: each function of the fundamental system (4.12) can be analytically continued to any point in C\{a 1 , . . . , a n } via any path that does not pass through any of the singularities. Concretely, if τ γ j denotes the operator of analytic continuation along the fundamental loop γ j , compare figure 2, then by problem 4.1 condition (b), any fundamental matrix Φ(z) satisfies But the same is also true for the matrix-valued function which, according to theorem 4.10, is invertible and bounded at z = ∞. Using Liouville's theorem and the piecewise constant, that is z-independent, form of G(z) it now follows that the matrix is single-valued on CP 1 and analytic everywhere expect at the singularities a 1 , . . . , a n ∈ C. Moreover, using problem 4.1, condition (c), (4.18) and our discussion on Vekua's regularization, it is clear that a 1 , . . . , a n are regular singular points of the system (4.19) in the following sense: Definition 4.11. Given a system (3.1) where A(z) ∈ C p×p is analytic in a punctured disk D r (z 0 )\{z 0 } for some r > 0, we say that z 0 is a regular singular point of (3.1) if any solution of the system has at most polynomial growth in a vicinity of z 0 .
Since the system (4.19) also has the given monodromy (G 1 , . . . , G n ), the above analysis solves RHP 2.1 in interpretation (b). Precisely Theorem 4.12  . Any matrix group with n ∈ Z 1 generators G 1 , . . . , G n ∈ GL(p, C) satisfying the constraint G 1 · . . . · G n = I can be realized as the monodromy group of a p × p linear system (3.1) on CP 1 having only regular singularities.
Plemelj's 1908 paper does not stop with theorem 4.12: indeed, while all singularities in a Fuchsian system are always regular singular points (this is Sauvage's 1886 theorem, see [27, theorem 2.1] or [83, theorem 16.10]), the converse statement is in general false. For instance, the non-Fuchsian system has a fundamental solution of the form 9 The choice of γ σ j intimately ties to the blow up constraint near z = a i in problem 4.1 (3), see [138, p 95]. and thus z = 0 is a regular singular point of system (4.20). Plemelj was aware of the distinction between regular singular and Fuchsian singular points, thus he applied a gauge procedure that took him from (4.19) to another system with equal monodromy and same singular points. Still, Plemelj's almost solution of RHP 3.1 introduced the idea of rephrasing a dynamical system, here (3.1) having only regular singularities, as a Hilbert boundary value problem. In turn, the analysis of the problem's underlying singular integral equations led to valuable information about the dynamical system itself, here theorem 4.12. We will encounter both themes throughout our discussions of OPSFA-type problems after the upcoming section.

Hilbert's 21st problem after 1908
Plemelj's 1908 result became widely accepted as providing a positive answer to RHP 3.1. Thus, in the following 75 years, the field connected to Hilbert's 21st problem shifted its focus towards the effective construction of Fuchsian systems with prescribed monodromy groups.
Here are the major results of this period: • In 1913, Birkhoff [14] revisited Plemelj's paper [121] and set out to simplify his argument based on successive approximations while also extending RHP 2.1 to certain difference equations. • In the late 1920s, Lappo-Danilevskiȋ [104,105] solved RHP 3.1 constructively, provided all generators G i are sufficiently close to the identity matrix. His method expressed solutions of a Fuchsian system and their associated monodromy via convergent series of the system's matrix coefficients. The solution of Hilbert's 21st problem subsequently boiled down to the problem of inverting the series and studying its convergence. • In 1956, Krylov [100] explicitly solved RHP 3.1 for all 2 × 2 systems with n = 3 singular points. His work made crucial use of Gauss hypergeometric functions. • In 1957, Röhrl [126] introduced a novel set of algebro-geometric ideas in the analysis of Hilbert's 21st problem. This allowed for a reformulation and generalization of RHP 2.1 to holomorphic vector bundles over Riemann surfaces, nowadays summarized under the umbrella of Riemann-Hilbert correspondences. Important early contributions to this field were achieved by Deligne [46], Kashiwara [91,92] and Mebkhout [110,111]. • In 1979, Dekkers [45] proved the solvability of RHP 3.1 in case p = 2 for an arbitrary number of singularities. His work did not directly address Hilbert's 21st problem (which was thought to have been solved by Plemelj), but contains its positive solution for p = 2 as special case. • In 1982, Erugin [55] considered RHP 3.1 for all 2 × 2 systems with n = 4 singular points. His work established a remarkable connection of this specialized problem and the Painlevé-VI equation. Compare [60, chapter 2, section 1.3] for more details and also [88].
Notice how all of the above works either yielded a positive solution to some special case of RHP 3.1 or were concerned with generalizations of the same problem to more general Riemann surfaces. However, as we already mentioned, RHP 3.1 is in general unsolvable as there are monodromy representations which cannot be representations of any Fuchsian systems. This surprising fact came to light in the important works of Bolibrukh in 1989, see [18,20]. In more detail, Bolibrukh's first counterexample concerns the following 3 × 3 system with n = 4 singular points: This system has Fuchsian singularities at a 1 = −1, a 2 = 1 2 , a 3 = 1 and a non-Fuchsian singularity at a 4 = 0. Although non-Fuchsian, Bolibrukh then showed that those singularities are regular singular points of (5.1) in the sense of definition 4.11. Moreover he proved that system (5.1) has non-trivial monodromy, but there exists no Fuchsian system with the same mondromy group and encoded singularity locations. Thus RHP 3.1 is in general unsolvable. The details of Bolibrukh's argument can be found in the monograph [5, chapter 2] and they are very different from Plemelj's analysis of the boundary value problem 4.1. Rather than discussing those techniques we only mention two features of Bolibrukh's counterexample (5.1). One, the sensitive dependence on the location of the singular points: once slightly perturbed, the answer to RHP 3.1 with the same monodromy can become positive, see [18]. Two, the mondromy representation of system (5.1) is reducible, a fact which intimately links to the following general result by Kostov [98] and Bolibrukh [19,20]: In the years following 1992 and theorem 5.1, Bolibrukh sharpened his results and derived a further series of sufficient conditions for the positive solvability of RHP 3.1, compare [21], these are conditions formulated in the language of holomorphic vector bundles and we refer the reader to the excellent monograph [5] for details.
In closing of this short section, and our content on the original RHP 2.1, we emphasize that Bolibrukh's application of methods from complex analytic geometry provided a definite, albeit negative, solution to Hilbert's 21st problem for Fuchsian systems on CP 1 . This problem was formulated back in 1900 and Bolibrukh lectured on his breakthroughs precisely 94 years after Hilbert at the ICM 1994 in Zürich, sadly a few years later he passed away.

Developments tangential to Plemelj's work-five examples
At this point we readjust our focus and turn away from Hilbert's 21st problem, the original RHP. Instead we highlight certain developments parallel to Plemelj's 1908 work on Hilbert boundary value problems and associated singular integral equations which found multiple applications in the OPSFA realm. Those appearances gave rise to an analytic apparatus, the Riemann-Hilbert techniques, which we are about to highlight through five examples-section 9 contains another example, but it is of different technical nature and more advanced. Overall, it is important to trace several of these techniques back to their origin [121], namely to Plemelj's work on RHP 2.1.

The Wiener-Hopf method
The first member of our analytical toolbox, a.k.a. the Riemann-Hilbert techniques, is the Wiener-Hopf method (pioneered by Wiener and Hopf in 1931 [140], but applied implicitly by Carleman before 1931). This method is commonly used in hydrodynamics, diffraction or linear elasticity theory and was originally developed for the investigation of the Milne-Schwarzschild integral equation which is used in models for radiative processes in astrophysics. Rather than discussing the Wiener-Hopf method in the more general setup we will simply use the original problem (6.1). How do we solve (6.1) with Riemann-Hilbert techniques? According to the Wiener-Hopf recipe one first extends (6.1) to the full real line via φ(x) := 0 for x 0 and obtains where ψ(x) := 0 for x > 0 and otherwise appropriately chosen for x 0 to achieve equality in (6.3). Then one takes a formal Fourier transform of (6.3) and finds The subscripts ± in (6.4) indicate our desire that Φ + (z), resp. Ψ − (z) admit analytic continuation to some part of the upper, resp. lower z-plane-in fact we will seek solutions of (6.1) which do not exponentially grow at x = +∞, and this implies that Φ + (z) is the limiting value of a function Φ(z) analytic for z > 0 and likewise Ψ − (z) the limiting value of a function Ψ(z) analytic for z < 0. Next, by elementary calculus, we find (6.5) and now Riemann-Hilbert techniques enter our calculation.
Hilbert Boundary Value Problem 6.1. Determine F(z) ∈ C such that (a) F(z) is analytic in C\R and extends continuously from either side to the real axis.
(b) On the real axis, the pointwise limits (c) Near z = ∞ we enforce the asymptotic behaviour as z → ±∞, the unique solution of problem 6.1 is readily seen to be by the Plemelj-Sokhotski formulae. The point of (6.6) is that we can use the underlying Wiener-Hopf factorization back in (6.4) and obtain Here, the equation's left-hand side defines a function which admits analytic continuation to the upper half-plane and the right-hand side a function which admits analytic continuation to the lower half-plane. Hence, there is an entire function E(z), z ∈ C such that as z → ∞, z > 0 and thus with problem 6.1 and Liouville's theorem, E(z) ≡ ic ∈ C. All together, This result at hand, we would now like to return to φ(x) via an inverse Fourier transform but how do we choose the contour Γ ⊂ C in (6.7) to achieve equality? Notice that by the Plemelj-Sokhotski formula, F + (0) = 1 √ 3 = 0, so the continuation of Φ(z) to the lower halfplane will have a double pole at z = 0 and a branch cut extending from −i to −i∞ along the imaginary axis, compare (6.5). Hence, if x < 0, we choose Γ as straight line in the upper half-plane and find that (6.7) vanishes by residue theorem. If x > 0, then we wrap Γ around the negative imaginary axis and encircle the double pole at z = 0 as well as the branch cut, compare figure 4.
Computing the residue and parametrizing the integral along the branch cut, one then obtains the well-known solution formula In summary of this short subsection, using scalar Hilbert boundary value problems we can solve integral equations of the form (6.2) by the Wiener-Hopf method. This method has found countless applications and is particularly well-suited for PDE problems with semi-infinite boundaries, e.g. the Sommerfeld diffraction problem, see for instance [114] and [33, chapter 16] for further details.

The integrable systems revolution of the late 1960s
The classical (think of Euler, Hamilton, Jacobi, Liouville, Neumann and Kowaleskaya) field of integrable systems found renewed interest in the late 1960s when it became clear that several nonlinear PDEs in 1 + 1 dimensions can be integrated by the inverse scattering method. Examples are the Korteweg-de Vries equation, the nonlinear Schrödinger equation (NLS), and the sine-Gordon equation, among many others. Important early contributions to this remarkable technique are due to Gardner et al [63], due to Lax [106], due to Zakharov and Faddeev [143], and due to Zakharov and Shabat [144]. For a solid first introduction to the inverse scattering method we refer the interested reader to the monographs [1,12,56,116], here we shall discuss only one example and showcase the appearance of Riemann-Hilbert techniques: consider the defocusing NLS How do we solve the corresponding initial value problem with Cauchy data y(x, 0) = y 0 (x) ∈ S(R) in the Schwartz space on the real line? Well, one first computes the reflection coefficient r(z) ∈ S(R) associated to y 0 through the direct scattering transform. There is not much freedom in this computation as the map y 0 → r is a bijection from S(R) onto S(R) ∩ {r : sup z∈R |r(z)| < 1}, see [11]. After that, one considers the so-called inverse scattering problem which is most conveniently formulated as the following Hilbert boundary value problem-this insight grew out of the works by Manakov, Shabat and Zakharov in the mid 1970s: is analytic in C\R and extends continuously from either side to the real axis.
(b) On the real axis, the pointwise limits where I ∈ GL(2, C) is the identity matrix.
Indeed, provided this problem is solvable, its (unique) solution solves the PDE (6.8) with initial condition y(x, 0) = y 0 (x) via the formula y(x, t) = 2iX 12 1 (x, t). (6.10) In order to arrive at (6.10) we use a slight modification of Plemelj's argument that lead him to the linear system (4.19): the jump matrix G(z), z ∈ Γ in problem 4.1 was piecewise constant whereas (6.9) presents us with the other extreme, namely a jump matrix which depends on z, x and t. However, assuming that problem 6.2 is solvable, we can define and then conclude that both, are single-valued, entire functions in z ∈ C. Moreover, using problem 6.2, condition (c) and Liouville's theorem we obtain the explicit formulae in term of the matrix coefficients X i j k = X i j k (x, t) and the third Pauli matrix σ 3 := 1 0 0 −1 .
The overdetermined system (6.11) forms one of the celebrated Lax pairs for the PDE (6.8), indeed writing out the Frobenius integrability condition for ∂Y ∂x = AY, ∂Y ∂t = BY, equivalently the compatibility condition we find entrywise the coupled partial differential equations However, the coefficients back in condition (c) of problem 6.2 satisfy trX 1 = 0 and Thus the above four coupled equations boil down to where we used (6.10) and so after substitution into one another which is the defocusing NLS (6.8). Summarizing this short computation, provided problem 6.2 is solvable, the linear Hilbert boundary value problem solves the nonlinear evolution equation (6.8) via the formula (6.10). This is a far reaching generalization of Plemelj's idea to construct solutions of certain linear ODE systems in terms of solutions of Hilbert boundary value problems, nevertheless the common approach to both problems is clearly visible. In section 7 below we will address the solvability question of problem 6.2 which requires a different toolset than the one used by Plemelj in his analysis of problem 4.1.

Painlevé special function theory
Painlevé functions form a family of special functions which is widely regarded as the substitute for classical special functions in nonlinear mathematical physics (such Airy, Bessel or Hypergeometric functions). Although Painlevé transcendents are non-expressible in terms of a finite number of contour integrals over elementary, elliptic or finite genus algebraic functions, several of their key analytic and asymptotic properties can be studied with Riemann-Hilbert techniques. These techniques are therefore in a way the analogues of contour integral representations and steepest descent asymptotic methods used in the analysis of classical special functions. For a comprehensive introduction to Painlevé special functions and Riemann-Hilbert techniques associated with them, including several references to their applications in mathematical physics, we recommend the two monographs [60,67]. Similar to the last subsection we will showcase the Riemann-Hilbert approach to one particular Painlevé equation, namely the homogeneous Painlevé-II equation, which, to a certain extent, is a nonlinear Airy equation. How do we solve the associated initial value problem? Well, in the Riemann-Hilbert approach to (6.12) one parametrizes the solutions u = u(x) of (6.12) not in terms of Cauchy data, but in terms of the monodromy data of an associated linear system of ordinary differential equations: precisely, define the generically two-dimensional real manifold and introduce with s k+3 = −s k , k = 1, 2, 3 for k = 1, . . . , 6 the triangular matrices Now consider the below Hilbert boundary value problem.

Hilbert Boundary Value Problem 6.3. For any x
(c) Near z = ∞ we enforce the asymptotic behaviour The route between problem 6.3 and (6.12) goes as follows: for any choice of parameters (s 1 , s 2 , s 3 ) ∈ S the Hilbert boundary value problem 6.3 is meromorphically with respect to x solvable, see [22]. In turn the unique solution to said problem leads to a solution of (6.12) via 12 1 (x, s), (6.13) which is the analogue of the NLS formula (6.10). Additionally, (6.13) satisfiesū(x) = u(x), i.e. the solution is real-valued on the real axis and conversely, every on the real axis real-valued solution of (6.12) admits a unique Riemann-Hilbert representation (6.13) for suitable s. In short, we do parametrize solutions of (6.12) in terms of the monodromy data (s 1 , s 2 , s 3 ) ∈ S using a bijection between the initial value solution space of (6.12) and S.
The derivation of (6.13), and thus the Riemann-Hilbert approach to the second Painlevé equation is due to Flaschka and Newell [57] and Jimbo et al [85][86][87] in the early 1980s. In a nutshell, we use again Plemelj's basic idea: suppose problem 6.3 is solvable and set

Then both,
are single-valued, entire functions in z ∈ C away from a discrete set D ⊂ R. Hence, using condition (c) in problem 6.3 and Liouville's theorem, we obtain the explicit formulae in terms of the matrix coefficients X k and from those, reading the compatibility condition entrywise, the following coupled ordinary differential equations, However problem 6.3 also possesses certain implicit symmetries that yield tr X 1 = 0, and which simplify the above four equations to where we used (6.13). Substituting both equations into one another we find and therefore Painlevé-II (6.12). In summary, the Hilbert boundary value problem 6.3 linearizes the nonlinear Painlevé-II ODE and therefore shares common ground with (4.19) and (6.8).
Before we move to our next example it will be important to mention that the monodromy data s = (s 1 , s 2 , s 3 ) ∈ S forms a complete set of first integrals for (6.12). Building on this feature, one can in principle use the highly transcendental equations in the further analysis of the Painlevé-II transcendents. This idea is at the heart of the isomonodromy method and the corresponding direct monodromy approach to Painlevé equations as pioneered by Its, Novokshenov, Kapaev and Kitaev in the mid 1980s, see [81] or [60, chapter 7]. Their approach relies on complex WKB techniques applied to problem 6.3 and certain a priori assumptions on the behaviour of u = u(x). We will later on highlight a different method, the inverse monodromy approach, a.k.a. Deift-Zhou method [37], which (in the extended form of Deift et al [44]) has become extremely popular in nonlinear mathematical physics since its discovery back in 1993. This method is suitable for the analysis of Painlevé functions and other problems that we are about to discuss.

The Heisenberg antiferromagnet
Our next example is concerned with a quantum mechanical model for an antiferromagnetic spin chain. It was introduced by Lieb et al [107] in 1961 and is known as the spin-1 2 XY model, equivalently the Heisenberg XX0 antiferromagnet. In this model 1 2 -spins are situated on a one-dimensional periodic and isotropic lattice, we allow only nearest neighbour interactions between the spins, neighbouring spins tend to point in opposite directions and a moderate external transverse magnetic field influences the spins. Subject to these four constraints the Hamiltonian of the model equals, see [107, (2.1)], where we set the nearest neighbour coupling constant to unity, 0 < h < 2 denotes the strength of the external magnetic field and are the Pauli spin matrices. The paper [107] computed the model's ground state, all elementary excitations and the free energy through creation-annihilation operator techniques. Here we shall focus on an important correlation function, the so-called emptiness formation probability P n := Prob {there are n ∈ Z 1 adjacent parallel spins up in the ground state} , which was evaluated in the thermodynamic limit N → ∞ for fixed h in [54, (10.2), (10. 3)] by the quantum inverse scattering method. The result is the Fredholm determinant formula where U n : L 2 (J, dx) → L 2 (J, dx) denotes the trace class operator with kernel acting on the interval J := (−Λ, Λ) ⊂ R with cosh 2Λ = 2 h > 1. The central analytical challenges associated with P n consist in accessing its large n asymptotic behaviour and in identifying an underlying integrable system for it. Both questions yield physically relevant information for the spin chain model, so how do we go about them? First, one factorizes the operator U n as where D n denotes multiplication (D n f )(z) := z − n 2 e λ(z) (z + i) with the invertible change of variables λ = λ(z) defined as E is the evaluation (E f )(z) := f(λ(z)) and K n : L 2 (Γ, |dz|) → L 2 (Γ, |dz|) the trace class operator with kernel acting on the arc Γ := z(J) ⊂ C defined in (6.17) and shown in figure 6. Second, one notices that (6.18) is a particular instance of an integrable kernel, i.e. a kernel of the form for some functions f j , g j ∈ L ∞ (Γ, |dz|), 1 j N. The associated integral operator K : L 2 (Γ, |dz|) → L 2 (Γ, |dz|) has many remarkable properties, most importantly (1 − K) −1 (if existent) can be computed in terms of the solution of a naturally associated Hilbert boundary value problem. This insight was formalized in the important 1990 paper [78] by Its et al, see [34] for a concise review, and reads in case of (6.18) as follows: Hilbert Boundary Value Problem 6.4. For any n ∈ Z 1 and 0 < h < 2, determine X(z) = X(z; n, h) ∈ C 2×2 such that (a) X(z) is analytic in C\Γ and extends continuously from either side to Γ\{ξ,ξ}. See figure 6 for the jump contour Γ and its endpoints ξ andξ. (b) On the arc Γ ⊂ C the pointwise limits (c) In a small neighbourhood of z = ξ and z =ξ we impose the blow-up constraint (d) Near z = ∞ we require the asymptotic normalization . Problem 6.4 and the emptiness formation probability (6.15) are related in the following way: as shown in [40, proposition 6.12], the Hilbert boundary value problem 6.4 is uniquely solvable for any n ∈ Z 1 , h ∈ (0, 2) and its solution leads to the below recursion for P n , P n+1 P n = X 11 (0; n, h), (6.19) in terms of the (11)-entry of the solution to problem 6.4. This identity is more complicated than the representation formulae (6.10), (6.13) and its derivation does not use a Lax pair as in the NLS or Painlevé-II example. Nevertheless, (6.19) enables us to use Riemann-Hilbert techniques in the analysis of the spin chain model, i.e. Plemelj's basic idea occurs once more.
Next we realize that where α n ⊗ β n is the finite rank integral operator on L 2 (Γ, |dz|) with kernel (α n ⊗ β n )(z, w) = α n (z)β n (w) in terms of α n (z) := 1 2πi z n 2 and β n (w) := w − n 2 −1 . But [78] showed that problem 6.4 is solvable if and only if (1 − K n ) −1 exists, so we can further simplify (6.15), and evaluate the remaining Fredholm determinant with the general theory of finite rank operators, At this point we use the general fact, see [78] or (7.5) below, that the solution of problem 6.4 can be expressed in integral form, similar to Plemelj's quasi-regular equation (4.4), with F(z) := (F 1 (z), F 2 (z)) , g(z) := (g 1 (z), g 2 (z)) and Setting z = 0 in the integral formula for X(z) and comparing with (6.20) yields P n+1 P n = X 11 (0; n, h), i.e. the representation (6.19). In summary, the Hilbert boundary value problem 6.4 linearizes to a certain extent (note that (6.19) is a formula for the discrete logarithmic derivative of P n and not P n itself ) the highly transcendental Fredholm determinant (6.15). We will display the usefulness of (6.19) after the next subsection.

Invariant random matrix models
Our fifth and last example in this section concerns the statistical behaviour of eigenvalues of a random matrix drawn from a particular ensemble. We will focus on an especially well-studied ensemble, the so called unitary ensemble M(n) of n × n Hermitian matrices equipped with the probability measure Here, dM is the Haar measure on M(n) R n 2 , N ∈ Z 1 a scaling parameter and V : R → R assumed to be real analytic satisfying the growth condition in order to ensure that (6.21) is a bona fide model. Measures of the form (6.21) appeared seemingly first in the work of Hurwitz in 1897 [73], they appeared in a paper by Wishart in 1928 [142] and in particular in Wigner's work in the 1950s [141] who used them to model the statistical properties of highly excited energy levels of complex nuclei, see [47] for an excellent historic review of the subject. A classical fact of the setup (6.21) is that the location of the eigenvalues of a matrix drawn from the unitary ensemble form a determinantal point process. This detail allows us to compute several of the key statistical properties of the random matrix model from the underlying Christoffel-Darboux kernel Indeed, focussing on one important statistical property of the model ( where K n,N is the finite rank operator with kernel (6.22). Formula (6.23) goes back to the seminal works of Dyson, Gaudin and Mehta in the 1970s and it is the analogue of the emptiness formation probability formula (6.15) in the spin chain model. Thus, in order to compute the likelihood of large eigenvalue gaps from (6.23), we need to carefully analyse the asymptotic behaviour of the kernel (6.22). However, in all but a few cases, the orthogonal polynomials that built the kernel are not known explicitly, so how can we effectively use (6.22)? Well, as it happens, the polynomials {p n (x)} ∞ n=0 and the kernel K n,N (x, y) itself can be characterized through a Hilbert boundary value problem. This remarkable fact was discovered by Fokas et al [58,59] in the early 1990s and the details are as follows: Hilbert Boundary Value Problem 6.5. For any n ∈ Z 0 and N ∈ Z 1 , determine X(z) = X(z; n, N) ∈ C 2×2 such that (c) Near z = ∞ we enforce the asymptotic behaviour Fokas, Its and Kitaev showed that this problem is uniquely solvable for a given n ∈ Z 0 if and only if the nth monic OP p n (x) exists. Moreover, in case when its solvability is ensured, then we have for any x, y ∈ R, p n (x) = X 11 (x; n, N) and 25) in terms of the (11)-entry of X(z) for the polynomial and the (21)-entry of the matrix product for the kernel. Here are the details: assuming the solvability of problem 6.5, the jump constraint (6.24) states that the first column of X(z) is an entire function, But the asymptotic normalization at z = ∞ enforces at the same time that and (6.27) Hence, X 11 (z) must be a monic polynomial of degree n and X 21 (z) a polynomial of degree at most n − 1 by Liouville's theorem. Returning with this information to (6.26) and using the Plemelj-Sokhotski formula we then find for z ∈ C\R, But geometric progression yields for z / ∈ R, and thus after comparison with (6.27) Since X 11 (z) has already been found to be a monic polynomial, it must therefore be the nth monic orthogonal polynomial p n (x) for the measure e −NV(x) dx on R, so the first identity in (6.25) follows. The kernel identity is slightly more complicated and involves X 21 (z): from the integral formula of X 22 (z) in (6.28), by geometric progression for z ∈ C\R, This shows that X 21 (z) must be proportional to p n−1 (z), i.e. we have for some γ n−1 ∈ C\{0}, and thus by orthogonality and the very definition of the sequence Next we note that any solution of problem 6.5 must be unimodular in the entire plane, i.e. detX(z) ≡ 1 for all z ∈ C. Indeed, the function detX(z) is analytic for z ∈ C\R and continuous on the closed upper and lower half-planes with det X + (z) = det X − (z), z ∈ R since det G(z) ≡ 1. Hence, det X(z) is entire and with det X(z) → 1 from problem 6.5, condition (c) we indeed arrive at det X(z) ≡ 1. This in mind, we now use the Christoffel-Darboux identity in (6.22), and our formulae X 11 (z) = p n (z), X 21 (z) = γ n−1 p n−1 (z) with γ n−1 h 2 n−1 = −2πi. A simple matrix multiplication and comparison with (6.29) leads to the second identity in (6.25). Once more we summarize our discussion: the Hilbert boundary value problem (6.5) allows us to access both, orthogonal polynomials and their Christoffel-Darboux kernel. In turn, the Riemann-Hilbert characterization allows us to rigorously analyse relevant statistical quantities in the random matrix model (6.21), provided we can efficiently derive a large n asymptotic expansion for X(z; n, N). The development of such an efficient scheme has been one of the many highlights over the past 30 years in the Riemann-Hilbert toolbox and we shall return to this important advancement after the next section.

Hilbert boundary value problems in L 2 -spaces
The boundary values problems encountered in the last section are different from Plemelj's initial problem 4.1 in that they include jump contours extending to infinity, jump contours with open ends and jump contours that self-intersect. It is therefore not obvious how Plemelj's solvability proof can be lifted to the problems of section 6. For this reason we now give a short overview of the relevant theory in which Plemelj's analysis of (4.4) in a space of continuous functions is extended to spaces of integrable functions-for a more thorough discussion we direct the interested reader to the articles [42,145] by Zhou and Deift, Zhou. Here are the basic two assumptions of this section: Assumption 7.1. Let Σ be a contour consisting of a finite union of smooth oriented curves in CP 1 with finitely many self-intersections.
Next, we derive the analogue of Plemelj's singular integral equation (4.4) for problem 7.3: if X ± (z) solve problem 7.3, then by the Plemelj-Sokhotski formula and condition (b), a.e.on Σ.
Having established in (7.4) the analogue of Plemelj's (4.4), how do we solve this equation? In short, we also rely on a Fredholm Alternative argument, similar to lemmas 4.6 and 4.7: let denote a pointwise factorization of G with functions (I ± V ± ) ±1 − I ∈ L 2 (Σ, |dz|) ∩ L ∞ (Σ, |dz|). This factorization allows us to rewrite the central integral equation (7.4) in the compact form a.e. on Σ : ρ(z) = I + (C w ρ) (z), with the Cauchy operators C ± of remark 7.4. In turn we have the following central result which states that the operator 1 − C w is injective if and only if a certain homogeneous version of problem 7.3 is only trivially solvable.

Proposition 7.6 ([145, proposition 3.3])
. ρ ∈ L 2 (Σ, |dz|) solves the homogeneous equation if and only if X ± (z) := ρ(z)(I ± V ± (z)), z ∈ Σ solves the homogeneous version of problem 7.3, i.e. the problem of determining X ± ∈ L 2 (Σ, |dz|) so that (a) There exists a C p×p -valued function H ∈ L 2 (Σ, |dz|) such that In order to state the Fredholm alternative needed in the solvability analysis of problem 7.3 we note that all our jump matrices encountered in this note possess much more regularity than what is assumed in assumption 7.2, in fact the matrices are piecewise smooth, admit local analytical continuations and at possible intersection points satisfy a cyclic constraint in the style of (3.4)-for instance in problem 6.3 the total monodromy at z = 0 is trivial 10 . These properties (see [60, pp 103-104] for a formalization) together with assumptions 7.1 and 7.2 ensure that the following far reaching generalizations of lemmas 4.6 and 4.7 hold true, at least for all matrix-valued Hilbert boundary value problems in this article.

Proposition 7.7 ([145, proposition 4.3]).
If the factors V ± (z) = V ± (z; ζ) in (7.6) depend on a parameter ζ analytically, then either (1 − C w ) −1 is meromorphic in ζ or 1 − C w is invertible for no ζ. Proposition 7.7 is an analytic Fredholm theorem and together with proposition 7.5 yields two central results, given that 1 − C w is a Fredholm operator of index zero, compare remark 4.4.

Theorem 7.9 (Zhou's vanishing lemma). The L 2 -Hilbert boundary value problem 7.3 is solvable if and only if the corresponding homogeneous version of the problem, see proposition 7.6, is only trivially solvable.
In the usual, smooth, setting of a Hilbert boundary value problem, the homogeneous version of it has the same analyticity and jump behaviour, but the asymptotic normalization at z = ∞ is replaced by This requirement is completely analogous to the behaviours which Plemelj enforced in his accompanying and associated problems, compare for instance problem 4.5.
Before moving to the next section, we will give one application of theorem 7.9 to problem 6.2 and en route verify that the defocusing NLS with Cauchy data y 0 (x) ∈ S(R) is solvable: suppose X(z) solves problem 6.2 but with condition (c) replaced by Define Y(z) := X(z)X † (z) for z ∈ C\R with X † denoting the conjugate transpose matrix of X. By property (a) in problem 6.2 we find that Y is analytic in the upper z-plane, continuous down to the real line and by (7.7) decays of O(z −2 ) as z → ∞ in the upper z-plane. Hence, by Cauchy's theorem, we have R Y + (z) dz = 0, and, adding to this its conjugate transpose, find in turn Now read off the diagonal entries in the last matrix equation and conclude that so by the fact that r ∈ S(R) ∩ {r : sup z∈R |r(z)| < 1}, see [11], we have X − (z) ≡ 0 on R by continuity of X − (z). In turn, from condition (b) in problem 6.2 also X + (z) ≡ 0 on R and thus all together X(z) is analytic in the upper z-plane, continuous in the closed upper z-plane and Using Carlson's theorem [130, theorem 5.1.2 and corollary 5.1.3], this allows us to conclude that X(z) ≡ 0 for Im z 0 and with the same logic for z 0 all together X(z) ≡ 0. By Zhou's vanishing lemma , i.e. theorem 7.9, we then deduce that problem 6.2 is solvable for any x ∈ R, t > 0 and so the solution of the initial value problem (6.8) with y 0 (x) ∈ S(R) exists. A similar application of theorem 7.9 to problem 6.3 does not yield the desired solvability result for general choices of s ∈ S. Indeed, for some choices s ∈ S, the solution X(z) = X(z; x, s) of problem 6.3 ceases to exist for certain x ∈ R and this is because the corresponding Painlevé-II transcendent u = u(x|s) has poles on the real axis, compare subsection 8.1. On the other hand, the solvability of problems 6.4 and 6.5 is always guaranteed and can be established without using theorem 7.9.

A Hilbert boundary value problem as Swiss army knife
By now we have seen several Hilbert boundary value problems in OPSFA related problems and we have acquainted ourselves to a certain extent with their abstract solvability theory in sections 4 and 7. Still, besides providing a positive solution to RHP 2.1 in interpretation (b), a reader unaccustomed to Hilbert boundary value problems might not yet grasp their usefulness: after all, although they seem to linearize nonlinear dynamical systems such as (6.8) and (6.12) or rephrase certain Fredholm determinants and orthogonal polynomials, those underlying problems (i.e. problems 6.2-6.5) are in all, but a few trivial, cases not explicitly solvable. So what is the whole point of a Hilbert boundary value problem? Well, the last 30 years have clearly shown that Hilbert boundary value problems really possess all the fundamental properties of a contour integral formula known, and appreciated, in classical special function theory. To the point, Hilbert boundary value problems underlie a large class of integrable models and as such allow us to • systematically derive dynamical systems (continuous ones, discrete ones or hybrids) for the quantities under consideration. This feature was first used by Plemelj in his work [121] on Hilbert's 21st problem. We have seen the same approach in action while discussing (6.8) and (6.12) and will further discuss it in the upcoming subsections. In these modern applications one must mention the pioneering roles of Shabat and Zakharov in the derivation of nonlinear partial differential equations from Hilbert boundary value problems, the role of Krein [99] in obtaining differential equations from integral equations and last, but not least, the role of Its [96, pp 377, 417] and his emphasis on studying correlation functions in statistical mechanics and quantum field theories via Hilbert boundary value problems and their associated differential equations. • analyse the models asymptotically in their thermodynamical limits. It is this feature which puts Hilbert boundary value problems on the same ground as their linear counterparts, i.e. contour integral formulae, and which has turned the Riemann-Hilbert approach to OPSFA problems into an unprecedented success story over the past 30 years. However the relevant asymptotic techniques did not grow over night: early progress-using Gelfand-Levitan type integral equations-on the asymptotic analysis of nonlinear wave equations that are solvable by the inverse scattering method was achieved in the early 1970s, namely byŠabat [128], Manakov [109] and Ablowitz and Newell [2]. These works were not always fully rigorous but the gaps were covered in the 1980s by Buslaev, Sukhanov, by Novokshenov and by Novokshenov, Sukhanov, see [39, p 182] for references. The first step of using directly a Hilbert boundary value problem for asymptotic questions seems to have originated in the works of Manakov [109] and Its [ However, aside from the NLS case, these works 'only' managed to asymptotically localize the initial Hilbert boundary value problems around certain special points. These points are the analogues of stationary phase points in the classical steepest descent method, but the model problems near them could not be explicitly solved, unfortunately. As it turned out, those local Hilbert boundary value problems are precisely the ones one faces in the isomonodromy deformation approach to the Painlevé transcendents as pioneered by Its, Novokshenov, Kapaev and Kitaev, see [81]. Thus, up to the early 1990s, one was able to use Hilbert boundary value problems in the asymptotic analysis of nonlinear wave and Painlevé equations, however one needed certain a priori information about the solution's behaviour. This pitfall was then bypassed in the groundbreaking method of Deift and Zhou [37] in 1993, the nonlinear steepest descent method. This method, in the 1997 extended version of Deift, Venakides, Zhou [44], has become the standard tool in the asymptotic analysis of Hilbert boundary value problems and all upcoming asymptotic results in this section have been originally derived from it or re-derived with it.
We will now summarize the technical essence of the Deift-Zhou nonlinear steepest descent method: in complete analogy to the classical steepest descent method, the nonlinear steepest descent method exploits the analytic and asymptotic properties of a Hilbert boundary value problem's jump matrix G(z). For instance, in case of problem 6.2 the jump matrix in (6.9) is highly oscillatory when x or t tend to infinity. Thus one would like to transform these oscillations to exponentially small contributions through the use of a contour deformation argument. And indeed, noticing that (6.9) admits the factorizations one may now deform the jump contour in problem 6.2 (according to the asymptotic regime at hand while taking into account the analytic properties of the reflection coefficient) and thus localize the problem near the relevant stationary point z 0 = − x 4t . Away from that point, the deformed Hilbert boundary value problem will be asymptotically well-behaved in the sense that its jump matrix converges to the identity in L 1 (Σ, |dz|) ∩ L 2 (Σ, |dz|) ∩ L ∞ (Σ, |dz|) as the parameters become large. We are thus facing a small norm problem away from the stationary point and such a problem is amenable to an iterative solution method. Here is the technical core of the argument: uniformly on compact subsets of C\Σ with C > 0.

Proof. Rewrite the integral equation (7.4) as
a.e. on Σ and abbreviate (recall remark 7.4) the operator However by assumption (8.1) and (7.2), in L 2 (Σ, |dz|) operator norm, which shows that that the integral equation is solvable in L 2 (Σ, |dz|) for sufficiently large t by the Neumann series. Moreover, since its right-hand side R(z; t) satisfies the unique solution ρ(·; t) − I ∈ L 2 (Σ, |dz|) of (7.4) also satisfies In summary, using proposition 7.5, the Hilbert boundary value problem 7.3 is uniquely solvable for sufficiently large t and we find from (7.5) in turn (this is the only time we use the L 1 (Σ, |dz|) bound in (8.1)), uniformly on compact subsets of C\Σ. This concludes the proof of theorem 8.1.
Near the stationary point(s) non-trivial technicalities of the nonlinear steepest descent method appear: one has to, either explicitly or implicitly, construct appropriate model functions in terms of classical or Painlevé special functions or, more generally, in terms of solutions to other dynamical systems that approximate the Hilbert boundary value problem asymptotically near the stationary point(s). Combining such local constructions with the aforementioned estimates away from the stationary point(s) one then typically arrives at a global small norm problem which is amenable to theorem 8.1 Remark 8.2. The above outlined Riemann-Hilbert nonlinear steepest descent techniques have been successfully applied in high performance numerical simulations, see the recent monograph [136] for a summary of research efforts in this direction.
Theory aside, we now discuss a series of impressive results which have been derived with Riemann-Hilbert techniques.

Connection formulae for Painlevé-II transcendents
Let us return to (6.12) and focus on the solution family parametrized by the monodromy data This one-parameter family is known as the Ablowitz-Segur family and it plays a dominant role in random matrix theory and integrable probability, see subsection 8.3 below and section 9.
Depending on the values of γ, the solution u = u(x|s) of (6.12) has different analytic and asymptotic properties, see figure 7. For one, the Ablowitz-Segur solution u = u(x|s) is smooth and bounded for x > 0 and if x → +∞, then which is proportional to the super exponential decay behaviour of the Airy function on the positive real axis (recall that (6.12) can be viewed as a nonlinear Airy equation). However, things change dramatically on the half ray x 0: if γ ∈ [0, 1), then u(x|s) is smooth, bounded and with the following oscillatory behaviour as x → −∞, Expansion (8.3) is due to Ablowitz and Segur [3,4] and was partially proven in [26,68] by Hastings, McLeod and Clarkson, using Gelfand-Levitan type integral equations. The first complete proof of the leading order in (8.3) is due to Its,Kapaev,Suleimanov and Kitaev [79,93,133] via the isomonodromy method. The first nonlinear steepest descent based proof of (8.3) is in the paper [38] by Deift and Zhou. Next, for γ = 1, the solution u(x|s) is still smooth for x 0 but now unbounded, in fact as which, to leading order, can be found in the work of Hastings and McLeod [68] and to all orders in [38], using once more the nonlinear steepest descent method. Finally, if γ > 1, the qualitative behaviour of u(x|s) changes completely: its smoothness is destroyed at finite x < 0 and solutions blow up, compare the below x → −∞ asymptotic expansion, which is uniform away from the zeros of the trigonometric function appearing in the denominator of (8.5). The singular expansion (8.5) is originally due to Kapaev [90] via isomonodromy techniques and was re-derived in [23] by nonlinear steepest descent techniques. Note that formulae (8.2)-(8.5) solve the connection problem for the underlying Painlevé-II transcendent, namely the problem of completely determining the asymptotic behaviour of u(x|s) as x → −∞, given the same description of u(x|s) as x → +∞, or vice versa. Such connection problems are standard exercises in classical special function theory and commonly solved based on contour integral formulae. However, the Ablowitz-Segur transcendent u(x|s) does not admit any such explicit formulae so it is a great achievement of the Riemann-Hilbert method to be able to solve this genuinely nonlinear connection problem in terms of classical special functions. Even better, one can uniformize the qualitatively different asymptotic behaviours (8.3)-(8.5) at x = −∞ through the use of Jacobi elliptic functions. This recent breakthrough is contained in the paper [24] and was made possible by nonlinear steepest descent techniques.

The Heisenberg antiferromagnet
We now return to (6.15) and discuss two central results for P n proven in [40, section 6]. First, the leading order asymptotic behaviour, ln P n = n 2 ln sin which holds uniformly on compact subsets of ( π 2 , π) φ. Hence, an overall Gaussian decay of P n with a decay rate that decreases with increasing external magnetic field. Second, the emptiness formation probability (6.15) really depends on two variables, n and h, equivalently n and ξ = e iφ , compare (6.17). As it turns out P n shares an intimate relation with the Painlevé-VI equation, more precisely P n is the Jimbo-Miwa-Ueno τ -function [85, section 5] of a particular Painlevé-VI equation when viewed in those variables. Here are the details, see [40, theorem 6.48]: set t :=ξ/ξ, then where y = y(t; n) satisfies the Painlevé-VI equation This impressive connection between P n and Painlevé SF theory puts the spin chain model firmly on integrable systems ground and as such should be compared with the famous Jimbo and Miwa 1980 result [84] on the diagonal Ising correlation functions and their relation to a σ-version of Painlevé-VI.
How does n behave statistically for large n? This question was raised by Ulam in 1961 [134] and he himself conjectured the first moment convergence which constituted a higher order simulation of E( n ) than (8.6). At this point of the story, in 1999, the Riemann-Hilbert/integrable systems community in the form of Baik et al [8] picked up the loose ends and proved the following spectacular result, settling en route (8.6), (8.7) and then some. 3 ([8, theorems 1.1 and 1.2]). Let u(x|s) denote the unique solution of (6.12) with the monodromy data s = (−i, 0, i), compare subsection 8.1. Define the distribution function Then, for any x ∈ R, and for any k ∈ Z 1 , The distribution function F(x) in (8.8), and it really is a distribution function, see our discussion on the analytic and asymptotic properties of u(x|s) in subsection 8.1, had appeared a few years prior to the Baik-Deift-Johansson theorem. Namely in the 1994 work of Tracy and Widom [135] where it was shown to capture the fluctuations of the largest eigenvalue of a matrix drawn from the ensemble (6.21) with quadratic external field, a.k.a. the Gaussian unitary ensemble (GUE), in the large n limit. The appearance of a random matrix theory distribution function in a problem of combinatorics was a big surprise at the time, however it turned out not to be an isolated curiosity. The Tracy-Widom distribution F(x) (8.8) (that name stuck after [8]) has become an ubiquitous limiting distribution that made its way to the financial markets, to wireless communication networks, stochastic growth processes and statistical mechanics, just to name a few areas. See [124] for an engaging outreach article on the Tracy-Widom law. In fact, twenty years after [8] it has become clear that the Tracy-Widom distribution is a modern day bell curve in that it describes the random fluctuations of a large class of integrable probabilistic models, see figure 8 for a graph of F(x) and F (x).
The proof of theorem 8.3 is beyond the material covered in this short note, it is a tour de force journey through (de)-Poissonization arguments, Toeplitz determinants and finally the asymptotic nonlinear steepest descent analysis of orthogonal polynomials on the unit circle. These polynomials can be characterized via a Hilbert boundary value problem in the style of problem 6.5, we refer the interested reader to the monographs [9,127].

Universality in invariant random matrix models
Our last discussion in this section concerns the global and local statistical eigenvalue behaviour in the random matrix model (6.21). First, through our assumptions placed on the external field V : R → R, the model's mean eigenvalue density 1 n K n,N (x, x), compare (6.22), has a limit whose support Σ V := {x ∈ R : ρ V (x) > 0} is a finite union of intervals. The global limiting behaviour (8.9) contains the most basic result in random matrix theory as special case, namely the Wigner semicircular law which asserts the almost sure convergence of the rescaled empirical spectral distribution in the GUE to the limiting distribution with density While the limiting density ρ V in (8.9) is determined by the potential V, the local eigenvalue statistics are to a large extent independent of V in the large n, N limit. Precisely, and this is the celebrated universality property, the local scaling behaviour of K n,N (x, y) is completely determined by the local characteristica of ρ V . For instance, in case of (6.21), only two qualitatively different scenarios arise generically, compare figure 9 for the GUE: if x * ∈ Σ V is such that ρ V (x * ) > 0, then uniformly on compact subsets of R λ, μ. This is the 1996 bulk universality result of Pastur and Shcherbina [118] which was proven with the resolvent method. Soon after that Bleher and Its [15] and Deift et al [43] fine tuned the Deift-Zhou nonlinear steepest descent method of section 8 to problem 6.5 and arrived at the same result. Especially [43] has become the standard reference for nonlinear steepest descent techniques used in the asymptotic analysis of orthogonal polynomials, and it has been extended in various directions. Those directions lead in turn to generalizations of the bulk universality class (8.10) and extensions to other types of orthogonal polynomial random matrix models, see section 10 below for a few references. The second case concerns the scaling limit of K n,n (x, y) at a boundary point x * ∈ ∂Σ V where ρ V typically vanishes square root like, certainly in the GUE. In this case (8.10) is replaced by the limit lim n→∞ 1 (cn) 2 3 K n,n x * + λ (cn) 2 3 , if x * is a right boundary point and with a sign change (λ, μ) → (−λ, −μ) in the left-hand side of (8.11) for a left boundary point. The limit is again uniform on compact subsets of R λ, μ, it is expressed in terms of the Airy function w = Ai(z) and c > 0 denotes some constant. The first proof of the soft edge universality (8.11) is due to Deift et al [43] and in more explicit form due to Deift and Gioev [36], in both cases based on nonlinear steepest descent techniques. In summary of this short subsection, the earliest universality results in random matrix theory were almost exclusively derived with Riemann-Hilbert nonlinear steepest descent techniques and thus paved the way for subsequent discoveries. Nowadays a wide array of different techniques is available in the derivation of RMT universality classes, some more suitable to certain ensembles than others. Most recently, the work of Erdős and Yau [53] put forth a toolbox useful in the universality analysis of Wigner random matrices, a broader class of symmetric matrices with independent entries. It is an open question whether the Wigner ensemble is accessible via Riemann-Hilbert techniques.

One more example
In this section we will venture into the world of operator-valued Hilbert boundary value problems, a topic which has received far less attention so far, especially from a mathematically rigorous perspective. The only exception is the work of Its and Kozlowski [80] and to a certain extent the paper [82] by Its and Slavnov which we will discuss towards the end. Here is our motivation to think about such problems: consider the stochastic heat equation (SHE) with distributional initial data Z(X, 0) = δ(X) and a Gaussian white noiseẆ =Ẇ (X, T).
Using the Hopf-Cole transformation H(X, T) := − ln Z(X, T) the SHE is mapped to the celebrated Kardar-Parisi-Zhang (KPZ) equation which has become the quintessential model for random surface growth processes with several remarkable connections to a number of different physical phenomena, see [28] for a detailed summary. One of the long-standing goals in the analysis of the KPZ equation was achieved in 2011 when Amir et al [6] derived the exact probability distribution of its solution with narrow wedge initial data. This allowed in turn to finally verify the celebrated KPZ scaling hypothesis for the KPZ equation itself. Here are the relevant details: 1 ([6, theorem 1.1]) . For any X ∈ R and T > 0 the Hopf-Cole solution H(X, T) = − ln Z(X, T) to the KPZ equation (9.1) with initial data Z(X, 0) = δ(X) has the probability distribution which is independent of X ∈ R and which equals In (9.2) we choose a Hankel contour Γ around R 0 , we abbreviate κ T := (T/2) 1/3 and K σ is the integral operator with kernel The explicit formula (9. Hence, the KPZ equation interpolates between two universality classes and (9.2) is therefore a crossover distribution. From a Riemann-Hilbert/integrable systems viewpoint we are interested in the Fredholm determinant entering the Amir-Corwin-Quastel formula (9.2). Integrating by parts we can rewrite its kernel in the form with the measure dσ(t) := σ (t)dt. Besides the particular choice of σ = σ T,z in (9.2), equality in (9.5) holds as soon as σ : R → C is smooth and Note that (9.5) includes the Fermi factor choice which plays an important role in the analysis of the Moshe-Neuberger-Shapiro random matrix ensemble [89, (14)], [108, (23)] and in the analysis of non-interacting quantum many-body systems [48, (158)]. Now given the formal-think of σ(t) = χ [0,∞) (t)-similarities between (9.5) and the standard Airy kernel (8.11) one might ask for an analogue of the Tracy-Widom Painlevé-II formula (8.8) for the Fredholm determinant of (9.5) on L 2 (s, ∞), perhaps even an underlying Lax pair which would further strengthen the connection between theorem 9.1 and integrable systems theory. The task of finding an analogue of (8.8) has been completed in [6, section 5] using the techniques of [135], see theorem 9.21 below. Here, we will employ Riemann-Hilbert methods and en route reproduce the Amir-Corwin-Quastel formula [6, proposition 5.2] while identifying an underlying Lax pair. Our approach is motivated by the duality between the kernels involved, namely i.e. the matrix-valued methods of subsection 6.4 will be lifted to operator-valued ones.

Technicalities up front
Throughout we are considering kernels K σ of the form (9.5) where dσ is a positive Borel probability measure on R for which all moments are finite, i.e.
and which is absolutely continuous with respect to the Lebesgue measure.
Definition 9.2. Let p ∈ Z 1 . We require the following abbreviations: (a) The direct sum Hilbert space equipped with the standard inner product f , g H p := p j=1 R f j (t)g j (t)dσ(t) and induced norm f H p := f , f H p . (b) The space L 2 (R, dσ; C p×p ) of p × p matrix-valued functions with entries in L 2 (R, dσ), equipped with the norm (c) The space I(H p ) of Hilbert-Schmidt integral operators on H p of the form with kernel K ∈ L 2 (R 2 , dσ ⊗ dσ; C p×p ). (d) The matrix identity operator I p on H p .
Next, fix some region Ω ⊂ C and consider K = K(z) ∈ I(H p ) which depends on an auxiliary variable z ∈ Ω. We now define the notion of an analytic integral operator, see [80, p 1781]. Definition 9.3. We say that K(z) ∈ I(H p ) with kernel K(z|x, y) is analytic in z ∈ Ω, if (a) for any (x, y) ∈ R 2 , the map z → K(z|x, y) ∈ C p×p is analytic in Ω. (b) for any z ∈ Ω, the map (x, y) → K(z|x, y) is in L 2 (R 2 , dσ ⊗ dσ; C p×p ).
Furthermore, if Σ ⊂ Ω ⊂ C is an oriented contour which consists of a finite union of smooth oriented curves in CP 1 with finitely many self-intersections, we will need boundary values of K(z) ∈ I(H p ), see [80, p 1782]. Definition 9.4. We say that an analytic in z ∈ Ω\Σ operator K(z) ∈ I(H p ) admits continuous boundary values K ± (z) ∈ I(H p ) on Σ ⊂ Σ with kernels K ± (z|x, y) if (a) for any (x, y) ∈ R 2 , the map z → K ± (z|x, y) ∈ C p×p is continuous on Σ . (b) for any (x, y) ∈ R 2 , the non-tangental limits lim λ→z K(λ|x, y) = K ± (z|x, y), λ ∈ ± side of Σ at z exist.
At this point we record two basic results about the kernel (9.5).
Proof. If σ(t) is the aforementioned Fermi factor, then the claim is proven in [89, proposition 1.1]. For general σ we argue as follows: write and obtain that for any x 1 , . . . , x n ∈ (s, ∞) and z 1 , . . . , z n ∈ C, n i, j=1 Thus, with K σ (x, y) = K σ (y, x), the function K σ (x, y) is Hermitian positive definite, so the claim follows provided we can show that see [129, theorem 2.12]. In order to get to (9.7) we compute the elementary integral, apply Fubini's theorem, and estimate for any y ∈ R, see [117, section 9.7], Hence (9.7) follows from our finite moment assumption placed on the Borel measure dσ.
Next to the trace class property of K σ we will also need the following result which is wellknown at least for the Airy operator induced by the kernel (8.11), i.e. for σ(t) = χ [0,∞) (t) in (9.5).
Proposition 9.6. For every s ∈ R, the self-adjoint operator K σ on L 2 (s, ∞) satisfies 0 Proof. We employ the standard Fourier transform trick, see [9, lemma 6.15], and therefore use First compute for any f ∈ L 2 (s, ∞), using the characteristic function of (s, ∞) ⊂ R. But from the Fourier representation, for any u ∈ R, and thus, since dσ is a probability measure, where we used Plancherel's theorem in the first and second equality. This shows all together that 0 K σ 1 and so, by self-adjointness, also K σ 1 in operator norm. Finally, since K σ is compact, it is sufficient to establish that 1 − K σ is injective in order to conclude that 1 − K σ is invertible on L 2 (s, ∞). So assume there is f ∈ L 2 (s, ∞)\{0} such that K σ f = f. Then all inequalities above must be equalities for this f, i.e. in particular But dσ is a positive Borel probability measure which is absolutely continuous with respect to the Lebesgue measure, so we must have σ − a.e. : and thus by the analytic properties of the Airy function, f = 0 ∈ L 2 (s, ∞), a contradiction.
We now combine lemma 9.5 with proposition 9.6 and conclude that the Fredholm determinant is well-defined on one hand and on the other hand satisfies 0 < F(s) < 1 for all s ∈ R. Even better, since whenever s < s , we also have that F σ (s) is strictly increasing in s, so F σ (s) is a distribution function of some random variable (one also needs Deift's lemma [6, lemma 2.20] for the regularity of F σ (s) in s). Additionally, see [131, theorem 3], there is a determinantal process with correlation kernel K σ . These features of (9.9) and (9.5) are not surprising in light of Johansson's discussion of the interpolating process in [89, section 1.4]. However we emphasize that [89] uses the Fermi factor choice for σ, here we are in a more general setup.

Solvability of the operator-valued problem
We now proceed to characterize (9.9) through a naturally associated Hilbert boundary value problem-very much as in subsection 6.4, but the problem is going to be operator-valued from the start. First, write where K σ,s has kernel K σ,s (x, y) := K σ (x + s, y + s). Then consider the following operators (we suppress the explicit s-dependence in our notation).
Definition 9.7. Let M i (z) ⊗ K j (z) ∈ I(H 1 ), i. j = 1, 2 denote the (0, ∞) z-parametric family of rank one integral operators with kernels and K j (z) are the integral operators on H 1 with kernel k j (z|y), Here is the relevant operator-valued Hilbert boundary value problem.
Our first result concerns the unique solvability of problem 9.8 Theorem 9.9. The Hilbert boundary value problem 9.8 is uniquely solvable for any s ∈ R and its solution takes the form (9.15) where N i (w) are the operators on H 1 which multiply by the functions n i (w|x) determined from the integral equation Proof. Suppose first that problem 9.8 is solvable. By condition (a), estimate (9.14) and the moment assumption on dσ, the Fredholm determinant is well-defined in the indicated z-domain by Hadamard's inequality (by condition (a) only, X 0 (z) is a Hilbert-Schmidt operator, but (d) ensures that the ordinary Fredholm determinant on H 2 converges). Moreover, by Morera's and Fubini's theorem, analyticity of X 0 (z) in C\[0, ∞) implies analyticity of e(z) in C\[0, ∞). Furthermore, using (9.14) again and the dominated convergence theorem, the boundary values e ± (z), z ∈ (0, ∞) satisfy e ± (z) = det X ± H 2 , z ∈ (0, ∞). (9.17) Our immediate goal is to derive a scalar Hilbert boundary value problem for e(z): from condi- we obtain from the Plemelj-Smithies formula [65, chapter II, theorem 3.1] that det(G(z) H 2 ) = 1 for all z ∈ R. With this in hand we return to (9.17), apply condition (b) and the multiplicativity property of Fredholm determinants, in turn Moreover, estimate (9.14) yields e(z) → 1 as z → ∞, so we are left to address the behaviour of e(z) in a vicinity of z = 0. By (9.12) and the estimate (9.13) we see that the operator P(z) := det I 2 + X reg (z) H 2 I 2 + X reg (z) −1 , |z| < on H 2 is well-defined even if det(I 2 + X reg (z) H 2 ) vanishes at some points in the disk |z| < , see [65, chapter VIII, theorem 1.1]. Hence from (9.12), as z → 0 and z / ∈ [0, ∞), ] . This shows that |e(z)| c|lnz| when z → 0 and so all together, e is analytic in C\{0} with a possible blow up at the origin that is logarithmic at worst. Hence, e(z) must be entire and with its limiting behaviour at z = ∞, thus e(z) ≡ 1 for all z ∈ C. As a consequence of this, we now know that any solution X(z) ∈ I(H 2 ) of problem 9.8 is invertible for all z ∈ C\[0, ∞) and so are its continuous boundary values X ± (z), z ∈ (0, ∞). Thus, if X (1) (z), X (2) (z) are two solutions of problem 9.8, then is analytic in C\[0, ∞), has continuous boundary values Y ± (z) on (0, ∞) and it obeys a problem in the style of problem 9.8 with Moreover, the kernel Y 0 (z|x, y) of Y 0 (z) has at worst a square logarithmic singularity at z = 0, so the map z → Y 0 (z|x, y) is entire (σ ⊗ σ)-almost everywhere. But Y 0 (z|x, y) → 0 as z → ∞ also (σ ⊗ σ)-almost everywhere, so by Liouville's theorem we conclude Y(z) = I 2 and which proves that problem 9.8, if solvable, is uniquely solvable. Second, we prove solvability of problem 9.8 and recall to this end that 1 − K σ,s is invertible on L 2 (0, ∞) for any s ∈ R, see proposition 9.6. Thus the integral equation (9.16) is uniquely solvable and by the boundedness of the resolvent we find for its solution, uniformly in x ∈ R, Now consider X(z) as in (9.15) and note that the kernel of X 0 (z) (the integral piece) is built out of the functions However, using (9.18) and standard estimates for Airy functions, see [117, section 9.7] we obtain by Cauchy-Schwarz Moreover, by the standard regularity properties of the resolvent operator, we also find from (9.16) that z → X i j 0 (z|x, y) is Hölder continuous on C\[0, ∞) for any (x, y) ∈ R 2 and i, j = 1, 2. Thus, by the Plemelj-Sokhostki formulae, z → X 0 (z|x, y) is analytic in C\[0, ∞) for any (x, y) ∈ R 2 and so the right-hand side in (9.15) analytic in z ∈ C\[0, ∞) according to definition 9.3. Having just established condition (a) in problem 9.8, we now turn to (b), (c) and (d): First, estimate (9.14) follows from (9.19) and the Plemelj-Sokhotski formulae. Second, z → X i j 0 (z|x, y) is Hölder continuous on [0, ∞) for any (x, y) ∈ R 2 , thus X(z) admits continuous boundary values X ± (z) ∈ I(H 2 ) on (0, ∞) according to definition 9.4 by the Plemelj-Sokhotski formulae. Now verify (9.11), so first compute from (9.15), On the other hand, by (9.12) and (9.15), and since, using here the general ( the chain of equalities (9.21) simplifies further But (9.16) is equivalent to the operator equation which, by symmetry of K σ,s , simplifies (9.22) to Condition (b) in problem 9.8 has therefore been established and we are left with the singular behaviour (9.12) and (9.13) near z = 0. Fix z / ∈ [0, ∞), 0 < < 1 and rewrite (9.15) as The remaining two integrals are operators in I(H 2 ) which are analytic in the disk |z| < and the same is true for the operator multiplied by the logarithm: indeed we already showed that and since the trace class G(z) ∈ I(H 2 ) is invertible with we also compute (as done previously for X − (z)G(z)), Hence, combining (9.24) and (9.25) with the explicit form (9.11), we deduce and from this, for z ∈ (0, ∞), the following operator composition formula for N i (z), (9.26) which is independent of the ± choice for the boundary values of X(z) on (0, ∞). Now use (9.26) back in (9.23) and conclude that N(z) near z = 0 is of the form i.e. since the logarithmic term drops out. All together, we can expand N i (z) near z = 0, which verifies the existence of the multiplication operators F i in (9.12) with their indicated properties, and estimate the kernel uniformly for w ∈ [0, ] ⊂ R, |z| < and any (x, y) ∈ R 2 . All combined in (9.23) we have now established (9.12) together with (9.13). This completes our proof.
During our upcoming computations we will use the following consequences of theorem 9.9.
Proof. The representation formula (9.33) is a general fact of integrable operators, see [34, (6), (7)]. In our case the finite sum in [34, (1)] is naturally replaced by a weighted integral but this has no significant impact on the formula for the resolvent kernel. In detail, since K σ,s is a Hilbert-Schmidt operator on L 2 (0, ∞), compare lemma 9.5, so is R σ,s = (1 − K σ,s ) −1 K σ,s and therefore R σ,s has a kernel R σ,s (x, y).
which yields for the underlying kernels, by (9.10), (9.16) and (9.28), The last corollary concludes our content on the operator-valued problem 9.8. We will now use this problem in the derivation of an integrable system for the Fredholm determinant F σ (s).

An operator-valued Lax pair
From now on it will be convenient, but not essential, to view the aforementioned multiplication operators M i (z), see definition 9.7, and N i (w), see theorem 9.9, as integral operators on H 1 with appropriate distributional kernels. To the point, we replace Subject to this convention we then compute, recall (9.26), which is the standard Airy system for the kernels and which translates into the integral operator identity Here, A i are z-independent integral operators on H 2 with distributional kernels to the explicit formula (9.15), write 1 w−z = − 1 z + w z(w−z) and obtain where X asy (z) ∈ I(H 2 ) and for any compact K ⊂ C with dist(0, K) > 0 there exists C > 0 such that for z ∈ C\K, X asy (z|x, y) C 1 + |z| 2 |s + x| 7 4 |s + y| 7 4 (9.38) uniformly in (x, y) ∈ R 2 . Similarly, from (9.27), and X asy (z) ∈ I(H 2 ) behaves similarly to (9.38). Inserting these two formulae for X(z) and (X(z)) −1 into (9.35) we obtain the exact operator identity with Y ∞ (z) ∈ I(H 2 ). Observe that for any compact K ⊂ C with dist(0, K) > 0 there exists C > 0 such that for z ∈ C\K, Y ∞ (z|x, y) C 1 + |z| |s + x| 7 4 |s + y| 7 4 , ∀ (x, y) ∈ R 2 . (9.39) Lemma 9.13. We have the finite rank integral operator identity valid for any i, j = 1, 2.
Proof. This follows from (9.30) by multiplying through with z and sending z → ∞ with z > 0, say. The double integral term vanishes in this limit by the dominated convergence theorem and the first term yields the stated identity by the same reasoning.
The symmetry lemma 9.13 allows us to simplify our expression for Y(z), we have now where A 0 is the integral operator on H 2 with kernel A 0 (x, y) = A i j 0 (x, y) 2 i, j=1 and The representation (9.40) mimics a Poincaré asymptotic expansion of Y(z) at z = ∞. We will now combine it with a corresponding representation at z = 0 which is computed from (9.23) and its counterpart for (9.27) when inserted into (9.35). The result equals with Y ori (z) ∈ I(H 2 ) analytic at z = 0 and where we have used the identity to ensure that all logarithmic terms cancel out in our computation of (9.41). At this point we combine proposition 9.12 with the representations (9.40), (9.41) and deduce from Liouville's theorem that y), (x, y) ∈ R 2 . Remark 9.14. The operator representation (9.42) is the analogue of [75, (9.4.28)] in the derivation of the Tracy-Widom Painlevé-II formula [135, (1.17)] through matrix-valued Riemann-Hilbert techniques. And, as in [75, (9.4.28)], the operator A −1 is nilpotent as a consequence of (9.33), Furthermore by (9.28), (9.16), for any z > 0, and which should be compared with the parametrization of the nilpotent matrix in [75, (9.4.30)].
Before continuing, let us summarize our computations as they will form the first half of the anticipated Lax pair. Proposition 9.15. Let N(z), z ∈ (0, ∞) denote the integral operator on H 2 as in (9.29). Then on H 2 have the kernels A i j k (x, y) stated above.
The second half of the Lax pair is computed in the same fashion, throughout replacing z-derivatives with s-derivatives. The main change occurs near z = 0 since the analogue of (9.35) with s-partial derivatives is analytic at z = 0. We thus have the simpler result Propositions 9.15 and 9.16 complete our derivation of the operator-valued Lax pair. We now proceed and show that this pair indeed produces an integrable system for the Fredholm determinant F σ (s), exactly the integro-differential Painlevé-II equation of Amir-Corwin-Quastel.

Derivation of the integro-differential equation
We first establish a connection between F σ (s) and problem 9.8.  (9.47) The coupled system (9.46), (9.47) is more complicated than its analogue [75, (9.4.34)] because of the aforementioned non-commutativity, still both systems share certain similarities. For instance so letting s → +∞ with the help of (9.10) and (9.16), using a Neumann series argument for n i (·|x) we find lim s→+∞ I(x, y) = 0 pointwise in (x, y) ∈ R 2 by the dominated convergence theorem. But I was already shown to be s-independent, thus it must be the zero operator as claimed.
We emphasize again that theorem 9.21 was first proven in [6] through an application of the Tracy-Widom method [135]. Here we obtained the same result from the operatorvalued Lax pair (9.43), (9.44) and the underlying Hilbert boundary value problem 9.8. Note that our method, en route, established the following novel features of u(s|t)-which should be compared with the ones for the Ablowitz-Segur transcendent u(x) = u(x|(−i, 0, i)) in subsection 8.1.
(a) the above boundary value problem for (9.55) is uniquely solvable, as u(s|t) = n 2 (0|t) is naturally constructed from (9.26), i.e. from the solution of the uniquely solvable Hilbert boundary value problem 9.8. (b) the map R s → u(s|t) is smooth for any t ∈ R in light of (9.49) and the analyticity properties of the resolvent. Equivalently because of (9.26) and the analyticity properties of X(z) = X(z; s), the solution of problem 9.8.

Remark 9.22.
A final remark is in order concerning the possibility to derive asymptotic expansions of F σ (s) as s → ±∞ directly from problem 9.8: while it is now fair to say that the operator-valued problem is well amenable to the rigorous derivation of dynamical systems, the development of a robust nonlinear steepest descent approach is work in progress by the author. As of now, the only rigorous operator-valued asymptotic techniques available are in [80,82], for two particular integral kernels. In that sense, operator-valued Hilbert boundary value problems are not yet on the same level as their matrix-valued counterparts. However, if one accepts other techniques, then certain asymptotic information for F σ (s) has been derived recently, albeit for very specific measures dσ, see [25,30,137].

Summary and further reading
This article is meant to provide a short survey of the theory of Riemann-Hilbert problems and their original appearances in mathematics. Particular emphasis is placed on Plemelj's 1908 work [121] on Hilbert's 21st problem and his use of Hilbert boundary value problems. Although Plemelj did not solve the problem for Fuchsian systems on CP 1 , his use of singular integral equations in the derivation of dynamical systems foreshadowed one of the key aspects of modern Riemann-Hilbert theory. Indeed, a Hilbert boundary value problem in modern terms is a Swiss army knife: once we know how a problem's observables can be characterized through a Hilbert boundary value problem then it presents us with a systematic tool to derive dynamical systems for the observables and with an efficient way to analyse them asymptotically-overall in complete analogy to classical contour integral methods. Still, the identification of a Riemann-Hilbert approach to a problem is non-systematic and we have displayed six different applications in the second half of the text. Much more than that has been done over the past decades and we would now like to give some guidance to the relevant Riemann-Hilbert literature-naturally based on the author's personal preferences: the monograph [60] is the standard reference concerning the Riemann-Hilbert approach to Painlevé equations and as a first taste to this subject we recommend the review articles [75,76] by Its. For a first introduction to OP Riemann-Hilbert methods in random matrix theory we single out Deift's notes [35], the notes by Bleher [16], by Its [77] and by Kuijlaars [101]. For continued reading on orthogonal polynomial methods in more general ensembles than (6.21) we recommend Kuijlaars' survey [102], the review article on universality [103] by the same author and the chapter on multi-matrix models by Bertola [13]. The standard reference to integrable operators is [34], the paper [40] and the somewhat non-rigorous monograph [96]. Regarding Ulam's problem and other combinatorial problems amenable to Riemann-Hilbert methods, the recent monograph [9] is likely going to become standard. Finally, and this is a topic which we did not touch upon, Hilbert boundary value problems have been identified as central tools in Gromov-Witten theories, i.e. in pure algebraic geometry. These theories aim to describe certain geometrical aspects of projective varieties in terms of invariants, suitably defined as intersection numbers on moduli spaces of curves. The invariants themselves are encoded in their generating functions whose convergence domains carry rich geometrical structures, nowadays known as Dubrovin-Frobenius manifolds. It was Dubrovin [49], see also [50], who first recognized the possibility to reconstruct such a manifold via a Hilbert boundary value problem. Further refinements of his procedure are given in the recent works of Cotti et al [31,32], sadly some of the last papers coauthored by Boris Dubrovin.