{'fraction_non_alphanumeric': 0.03880605640314378, 'fraction_numerical': 0.011760286638927416, 'mean_word_length': 4.220060331825038, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 13, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
QCD Jet Rates with the Inclusive Generalized kt Algorithms

We derive generating functions, valid to next-to-double logarithmic accuracy, for QCD jet rates according to the inclusive forms of the kt, Cambridge/Aachen and anti-kt algorithms, which are equivalent at this level of accuracy. We compare the analytical results with jet rates and average jet multiplicities from the SHERPA event generator, and study the transition between Poisson-like and staircase-like behaviour of jet ratios.


Introduction
The production of hadronic jets in high-energy collisions is one of the most striking and universal features of particle physics. It is now well understood that jets are associated with hard parton emission, followed by parton showering and hadronization. However, the observed features of jets depend not only on their intrinsic properties but also on the algorithms used to find them. These algorithms have been the object of many years of development and refinement -see [1] for a review.
The jet algorithms used most widely nowadays are of the sequential recombination type. Of these, the k t algorithm [2][3][4] was used extensively at LEP and HERA, and to a limited extent at the Tevatron. The good theoretical properties of the k t algorithm mean that jet rates and multiplicities can be calculated perturbatively, either in fixed order or in an all-orders logarithmic approximation. However, the irregular angular shapes of k t -jets make them less suitable for analyses of hadron collider events, where multiple parton interactions give rise to underlying activity not associated with the primary hard process.
A variant of the k t algorithm that allows better control of the underlying event is anti-k t [5], which has become the preferred jet finder at the LHC. It belongs to the general category of inclusive generalized k t jet algorithms, defined in [6]. Many of the good theoretical features of k t extend more broadly to the members of this category. In particular their leading and next-to-leading double logarithmic jet rates can be calculated to all orders, and are in fact the same for all members, including k t and anti-k t .
In the present paper we calculate these jet rates, along with the average jet multiplicities, in a variety of approximations. We do so by finding, in each case, the equations that the quark and gluon jet generating functions [7][8][9] must satisfy. In some cases, an explicit solution in terms of special functions can be found, while in others we are forced to resort to numerical methods (though we can get some analytic understanding of the properties of the solution, as we describe in detail for one case in an Appendix). We then compare the results of our theoretical calculations with those obtained via Monte Carlo simulations using the SHERPA event generator [10,11]. Our theoretical calculations of jet rates are carried out for the case of e + e − collisions, but we also perform calculations and simulations of the average jet multiplicity in hadron-hadron collisions, finding that similar behaviour is obtained.
The outline is as follows. We begin by recalling the definitions of the inclusive jet algorithms and the relevant Sudakov form factors. In Section 4, we introduce the jet generating functions and derive the integral equations that they satisfy at double-leading-logarithmic accuracy (DLA) and next-to-double-leading-logarithmic accuracy (NDLA). In Sections 5 and 6, we compute the resulting jet rates and average jet multiplicities, both at fixed coupling and also at next-to-leading order in the running coupling. In Sections 7 and 8, we compare with results of simulations performed using SHERPA for e + e − and pp collisions, respectively. In Section 9 we assess the ability for our analytic results to describe sub-jet multiplicities in boosted events. In Section 10 we consider the implications for the scaling patterns of jet multiplicities. Our conclusions are summarized in Section 11. Details of the derivation and properties of the partial differential equation (PDE) for the average jet multiplicity are relegated to appendices.

The inclusive generalized k t jet algorithms
We consider first the case of multijet production in e + e − annihilation, for which the inclusive algorithms are defined as described in the FastJet user manual [6], Sect. 4.5. The distance measures are with p = 1, 0, −1 for the k t [2], Cambridge/Aachen [12,13] and anti-k t [5] algorithms, respectively. At any stage of clustering, if a d ij is the smallest measure we combine objects i and j. If d iB is the smallest we call i a jet candidate and remove it from the clustering list. We then call jet candidates with energy E i > E R resolved jets. Thus the jet rates, at a given value of the e + e − centre-of-mass energy E cm , are functions of the radius parameter R and the minimum jet energy E R . This is in contrast to the exclusive k t (Durham) algorithm, where one effectively sets ξ R = 1 2 and continues clustering objects until all d ij are above some fixed value d cut = y cut E 2 cm , so that the jet rates are functions of the single parameter y cut .
In hadron-hadron collisions, the c.m. frame of the parton-parton hard scattering process is not known and therefore one has to adopt a longitudinally invariant form of the algorithms [3,4]. To that end, Eqs. (2.1) are replaced by where p ti , y i and φ i are the transverse momentum, rapidity and azimuth of object i, respectively, and we define jet candidates with p ti > E R as resolved jets. As far as leading logarithms are concerned, the jet rates defined by (2.1) and (2.2) will be the same, and therefore in the following Sections we refer mainly to e + e − annihilation. By "leading logarithms" here we always mean leading double and next-to-double logarithms, α n S log 2n and α n S log 2n−1 , where the logarithms are those of 1/R and/or Q/E R , Q being the hard process scale. By taking E R sufficiently large in hadron-hadron collisions, we avoid such leading contributions from initial-state showering and the underlying event, so these terms are determined by the timelike showering of final-state partons.

Sudakov factors
The evolution scale for coherent parton showering is ξ ≡ 1−cos θ with θ the emission angle. To be resolved, an emission must have ξ > ξ R and E > E R . The probability for a single resolvable gluon emission from a quark of energy E at scale ξ is thus where the running coupling is evaluated at the transverse momentum scale of the emission, in terms of the coupling at the hard scale, we have to next-to-double-log accuracy (NDLA) and Then the probability for no resolvable emissions (the quark Sudakov factor) is Similarly for a gluon, the probability of a single resolvable gluon, quark or antiquark emission is which gives to NDLA and the gluon Sudakov factor is Note that all this is independent of the value of p, so that all the inclusive generalized k t algorithms are equivalent at this level of precision.

Generating functions
By definition the generating function for resolved jets from a quark (i = q) or gluon where R i n is the corresponding n-jet rate, i.e. the probability of finding n resolved jets. Thus the jet rates can be recovered from the generating function by successive differentiation at u = 0: On the other hand the average multiplicity of resolved jets is obtained by differentiating at u = 1. Writing the average jet multiplicity from a quark or gluon of energy E at scale ξ as N i (E, ξ), we have The generating functions Φ q,g must thus satify the boundary condition The generating function for e + e − annihilation at c.m. energy E cm is that for two quarks of energy E cm /2, each filling one hemisphere: (4.5)

Next to double leading logarithms
For ξ > ξ R and E > E R , we have to NDLA Using the expressions given earlier to eliminate the Sudakov factors, we may write these in the equivalent forms The solution for the quark generating function is then easily seen to be We can solve for the gluon generating function by iteration, and then substitute in this equation to get the complete solution.

Leading double logarithms
In the DLA we keep only the singular parts of P gq and P gg and can drop P qg . For brevity we define the logarithms as Then ∆ q,g (κ, λ) = e −aq,gκλ (4.10) where a q,g = C F,A α S , and We can simplify the equation for Φ g by noting that where the boundary condition Φ g (u, κ, 0) = u for all κ ≥ 0 implies C(u, κ) = ln u. Thus (4.14) Comparing with (4.11) we see that in the DLA

Jet rates
From the quark and gluon generating functions in the previous section it is relatively straightforward to construct the jet rates to next-to-double log accuracy (NDLA).
In terms of the logarithmic variables (4.9) the quark and gluon integrated splitting kernels (3.4) and (3.7) are simply

(5.4)
In the fixed-order expansion it is not necessary to include the running coupling for the finite terms in the splitting functions, as doing so would affect results beyond the NDLA. In these finite term we can set terms proportional to e −κ to 0 after the integrations, which is equivalent to allowing the original z integration to range over (0, 1), since it is not singular in energy. However, we have to bear in mind that the relevant integrals vanish when κ = 0.
we solve the gluon generating function by iteration to third order in u.
In the result for Φ g we have definedκ = κ(κ ) when the κ integral ranges from 0 to κ (κ ). Defining the generating function for resummed rates up to 5 jets at NDLA is then Substituting in (5.4) and using (4.5), we find the rates for e + e − annihilation in the form for n ≤ 5 and j ≤ 3, where n is the number of resolved jets and R n,i has i powers of the logarithms (either κ or λ): One check on these results is that the DLA coefficients agree with the previous computation in [14] for R 44 and R 56 . A second is that the 2-jet inclusive fraction obtained by summing the rates is 1.

Average jet multiplicity
The average jet multiplicity, i.e. the mean number of resolved jets, as a function of the hard process scale Q, the angular resolution R and the minimum jet energy E R , provides a useful overall measure of jet activity and substructure. As indicated by Eq. (4.3), this quantity is obtained simply from the first derivative of the relevant generating function.

Next to double leading logarithms
Writing the average jet multiplicity in e + e − annihilation as N ee , from the fixed-order jet rates (5.13)-(5.16) we obtain to O(α 3 S ) Here the terms in b 0 originate from including the running coupling. We see that these terms enhance the average jet multiplicity with respect to a fixed-coupling calculation.
To perform an all-orders resummation to NDLA, we repeat the analysis of [15], for the inclusive algorithms instead of the exclusive k t algorithm. In terms of the generating functions, we have where N q,g , the average quark and gluon jet multiplicities, satisfy the equations In Appendix A we show that (6.3) is equivalent to the following PDE in terms of the logarithmic variables (4.9): Similarly we find from (6.4) where to the required accuracy we may set in the last term so that finally Note that the n f dependence in (6.8) is very weak and vanishes in the large-N limit: This is because at large N a qq pair from gluon splitting radiates like a gluon.

Leading double logarithms
Dropping the non-singular parts of the splitting functions, we have from (4.15) In the leading double log approximation, α S is treated as a constant. Then the solution to (6.11) is a modified Bessel function: The asymptotic behaviour for large argument, implies that for high energy and small cone size (6.14)

Running coupling
Taking into account the running of α S to next-to-leading order, we have Thus if we drop terms of relative order α 2 S the solution to (6.11) is which agrees with the b 0 -dependent terms in (6.1). However, for large κ and/or λ, b 0 (2κ + λ)α S ∼ 1 and therefore we need to take into account the running of α S to all orders.

Numerical solution: DLA
Treating the running of α S to all orders, as in Eq. (5.3), but still neglecting the finite parts of the splitting functions, we have in place of (6.11) This PDE is not straightforward to solve explicitly. Its properties are discussed in Appendix B. It may however be solved numerically by discretization. Writing and Equating these expressions, one can solve iteratively for g m+1,n+1 starting from the boundary values g 0,n = g m,0 = 1.

Numerical solution: NDLA
To include the finite parts of the splitting functions, we may write (6.5) and (6.8) with equivalent precision as The PDEs (6.21) can be solved numerically by a simple extension of the method outlined above. For the discretized κ-derivative, we use We can then write the right-hand side of (6.21) as and equate this to (6.19). Similarly, to obtain the quark jet multiplicity we write equate this to (6.24) with c g , δ g replaced by to obtain the discrete equivalent of (6.5), and solve iteratively for f m+1,n+1 starting from the boundary values f 0,n = f m,0 = 1.

Monte
Carlo results: e + e − In Fig. 1 results for the average jet multiplicity from the SHERPA Monte Carlo are compared with the fixed-coupling DLA results (for fixed α S = 0.118) and the NDLA results for running α S , the latter derived numerically as explained in the previous Section. The average jet multiplicity is shown for E cm = 200 TeV and E R = 10 GeV as a function of the jet radius parameter R. For these extreme values we expect the leading logarithms to dominate. At such a high energy, the asymptotic approximation (6.14) (green dashed) agrees well with the exact DLA formula (6.12) (solid green), even at R ∼ 1.
A small QCD scale Λ ∼ 50 MeV in the NDLA formulae gives reasonable agreement with the Monte Carlo results, in which Λ M S = 180 MeV. Such a scale change is a next-to-NDLA effect and therefore within the uncertainties. The Monte Carlo results follow the formulae up to ln(1/R) ∼ 4, after which the parton-level Monte Carlo result approaches the parton multiplicity (about 150 at the shower cutoff Q 0 ∼ 1 GeV), while the NDLA diverges towards the Landau pole of the running coupling at ln(1/R) ∼ 5. Surprisingly, the hadron-level result follows the fixed-coupling DLA perturbative prediction further, up to around 100 jets at ln(1/R) ∼ 5, before breaking away to saturate at the hadron multiplicity. We note that the quantitative agreement between Monte Carlo simulation, implementing the running coupling, and the analytic DLA result depends on the choice for the fixed coupling in the latter. For the extreme c.m. energy we consider, a somewhat smaller value for α S can improve the agreement with the simulation for larger values of ln(1/R). At the same time for lower values of ln(1/R), where we do not expect the DLA to fully apply, the agreement would get worse. Figures 2 and 3 show SHERPA parton and hadron level results for the jet multiplicity as a function of c.m. energy, compared with the DLA and NDLA predictions. Again, a QCD scale Λ ∼ 50 MeV in the NDLA formulae gives reasonable agreement with the Monte Carlo results. At the smallest values of R, the Monte Carlo values are approaching the NDLA predictions from below, indicating further subleading effects. As in Fig. 1, the agreement at very small R is better at hadron level than at parton level.
It can be see from both figures that the Monte Carlo results for all three inclusive algorithms, anti-k t , C/A and k t , are very similar, consistent with the lack of dependence of our analytical predictions on the power p in Eqs. (2.1) and (2.2). The jet multiplicity is systematically slightly lower for the k t algorithm, which we conjecture is due to its tendency to gather more low-momentum particles into jets from angles somewhat larger than the canonical jet radius R. This leads to a smaller number of jets for a given final-state multiplicity. This higher susceptibility of k t jets to radiation was discussed in [16] where also a first-order estimate for the resulting effective area of k t and C/A jets was derived. Qualitatively we observe the same behaviour as for the passive jet area discussed there: the effect increases for higher jet energies and larger nominal R.

Monte Carlo results: pp
The motivation for studying the inclusive family of jet algorithms is the similarity to phenomenologically relevant hadron collider algorithms. For jet rates, and therefore also the mean number of jets, the factorization of the PDF in the initial state holds at the double leading logarithmic order. In this section we compare the predictions of our generating functions with Monte Carlo calculations and determine the reliability of our resummation in the LHC context. Figure 4 shows results for pp collisions at the current LHC energy of 8 TeV. Here the average multiplicity of jets with p T > E R = 20 GeV is plotted as a function of H T /2, where the scalar sum of all jet transverse momenta. This choice corresponds closely with the initial parton transverse momenta in the QCD 2 → 2 hard scattering. Each point in the analytic result is weighted by the partonic fraction of final-state quarks versus gluons. In other words, the analytic result is where c q + c g = 1 and are determined from the proportions of final-state jets in the contributing 2 → 2 hard subprocesses. The considered values for R are the same as in Fig. 2. The multiplicity levels off and even decreases at high values of the primary parton p T , owing to the transition from gluon to quark jets. Overall we observe a good agreement of the analytic estimate with the SHERPA simulation. As before we use a fixed coupling of α S = 0.118 in the DLA calculation. In particular for large values of H T /2 a reduced value could be more appropriate, improving the agreement with simulation. Notably, the dependence on the jet algorithm for the simulated results is rather mild. As observed for e + e − collisions in the above, the k t algorithm tends to produce slightly less jets for a given final-state multiplicity.

Sub-jet multiplicities
Jet substructure is an increasing relevant handle for distinguishing strongly decaying new physics from the QCD background [17]. A particular observable of interest is the sub-jet (R (sub-jet) ∼ 0.1 − 0.4) multiplicity inside a larger fat jet (typically with a R (fat) ∼ 1−1.2) [18,19]. In this section we briefly outline how our generating function can predict QCD sub-jet multiplicities, and how these vary significantly from new physics signals.
In the first case, we modify our formula for the average jet multiplicity to account for the reduced available phase space of the fat-jet. At the DLA, the phase space boundary in the ξ integration is now the fat-jet radius R (fat) as opposed to the entire hemisphere, so that the angular logarithm in our average jet multiplicity (6.12)  becomes log(R (fat) /R (sub-jet) ). Further, we choose the fat-jet p T as the upper energy scale in analytic formula. With these modifications we compare the analytic formula with a SHERPA sample of pure quark jets in Fig. 5. The DLA result captures the scaling with respect to fat-jet p T , R (fat) and R (sub-jet) quite well. Even so we have to note, that the level of agreement is highly dependent on the chosen value for α S . As a second example we contrast in Fig. 6 our DLA multiplicity estimate with the result obtained for the production of a pair of hadronically decaying W bosons of mass 5 TeV. We again consider a centre-of-mass energy of 200 TeV resulting in rather boosted, collimated decays for the bosons. As a consequence we observe that in particular for sizable values of R, where the substructure of the W decays is not yet resolved, we expect many fewer jets than predicted for pure QCD production. As we decrease R we resolve more and more jets again, the hadron level jet-multiplicity being significantly higher than the one at parton level.

Jet scaling patterns
The work in [20][21][22] discussed the potential for extrapolating jet rates based on universal scaling patterns. These patterns are most easily classified in terms of the ratio of exclusive jet rates R (n+1)/n = σ n+1 /σ n . In the Durham (exclusive k t ) algorithm, it was found that for low multiplicity n ≤ n jets , emissions are essentially Poisson-like so that R (n+1)/n ∼ (n + 1) −1 . The tail of the multiplicity distribution then produces dominantly staircase or geometric scaling where R (n+1)/n ∼ constant. This regime is driven by the fractal nature of QCD radiation in the gluon dominated limit. We know from previous work that the expected scaling patterns of jets can depend dramatically on the jet algorithm. One example of this is the JADE algorithm, where the non-exponentiation of the primary emissions precludes the Poisson extrapolation even in the pseudo-abelian limit [23]. In this section we would like to address scaling in the inclusive generalized k t class of algorithms. This extends the results in [20] and strengthens the case for investigation at hadron colliders.

Poisson breaking components
With the leading logarithmic coefficients from Eqs. (5.13)-(5.16) it is easy to make some first statements about scaling in the generalized algorithm. It is clear from the structure of the coefficients that in the limit C A → 0 a perfect Poisson distribution emerges. Now a simple comparison between the generalized and Durham algorithms is the relative size of the Poisson breaking components in the lower multiplicity rates, for example the 2-gluon correlated emission contribution to the 4-jet rate (5.15). For the double-leading logarithmic coefficients to the 4-jet rate, R 44 , we find C Durham Hadron Level Parton Level normalization of this work. The lowest-order Poisson breaking term is relatively larger in the generalized class of algorithms. We would thus expect the onset of staircase (geometric) scaling to come about for even smaller values of the logarithm, and to better match the staircase behaviour at lower multiplicity. Evidence of this may be found in the Z + jets analysis presented in [20]. We present in table 1 the relative sizes of the Poisson breaking terms in the DLA expanding coefficients for the two algorithms compared with the idealized Poisson and staircase predictions. In this table we compute the 6-jet rate in the generalized algorithm by expanding the resummed jet rates from Eq. (5.11) to O(α 4 s ) which entirely determines the DLA resolved component.

R-dependence of scaling
A second question we wish to answer in this section is how the idealized scaling patterns depend on the jet radius parameter R in the generalized k t algorithm. The additional handle provided by the separation of the angular and energy regulation  allows us to probe an effect unbeknownst in the Durham algorithm. At the double leading logarithmic level, the resummed rates in the general algorithm are invariant under exchange of κ ↔ λ. Decreasing the size of the jet merely increases the overall logarithm.
In the simulation however, although we increase the overall logarithm when we require smaller resolved jets, we also change the relative contributions of the primary and secondary contributions due to kinematics. This effect is present even at the level of the ordered 2-gluon emission matrix element [24].
Using the parton shower we find that larger jet sizes dramatically increases the goodness of the fit with respect to the Poisson hypothesis. This has a clear interpretation. Correlated emissions in these events are predominantly intra-jet evolution for large jet radius. In addition, this effect is larger than that due to the breaking of κ ↔ λ symmetry at NDLA.
To analyze this further we investigate the observable θ 34 , defined as cos θ 34 ≡ p 3 · p 4 |p 3 ||p 4 | , (10.1) in reconstructed 4-jet events. The third and fourth hardest jets overwhelmingly originate from emitted gluons. As we suspect that in the simulation secondary Poisson breaking emissions are on average closer in θ 34 , this observable gives another estimate on the size of these breaking effects. Comparing the prediction of pseudoabelian QCD (in practice the parton shower with the g → gg and g → qq branching off) we see in Fig

Conclusions
In the present paper we have derived the generating functions of jet rates for the inclusive generalized k t jet algorithms, valid in the next-to-double log approximation (NDLA), and used them to compute jet rates and average jet multiplicities as functions of the jet radius parameter R and the minimum jet energy E R . At this level of precision, the results are independent of the power p that distinguishes the inclusive k t , Cambridge/Aachen and anti-k t algorithms.
The analytical results on e + e − annihilation are in broad agreement with those from the SHERPA Monte Carlo event generator, including the weak p-dependence. Surprisingly, the hadron-level Monte Carlo results follow the analytical predictions down to very small jet radii, well beyond the range of the perturbative parton shower, indicating that the cluster hadronization model in SHERPA is smoothly matched to perturbation theory. Analytical predictions for pp collisions at the current LHC energy also agree fairly well with SHERPA Monte Carlo results.
The SHERPA Monte Carlo results and generating function in the generalized k t algorithm were also used to study the transition from Poisson-like to staircase-like behaviour in the (n + 1)/n jet ratios. The relatively larger non-abelian terms in the DLA e + e − jet rates, compared to the Durham algorithm, indicate a transition at lower values of n. Furthermore, the jet radius was found to have a significant impact on scaling, where smaller jet sizes receive larger secondary contributions. This feature may be relevant for jet activity studies in jet substructure analyses.
These results represent the most detailed comparison to date between analytical and Monte Carlo predictions for the inclusive jet algorithms in current use at the LHC. Information about the dependence of jet rates and multiplicities at high jet energies and relatively small radii are particularly important for studies of boosted jets and jet substructure, which play an increasingly important role in searches for new physics at the LHC. As we demonstrated, our leading-logarithmic predictions are relevant to such studies, and their fairly good agreement with Monte Carlo results is encouraging. Possible improvements such as optimized scale choices and empirical subleading contributions would be worth exploring. Studies of other features of the jet multiplicity distribution along these lines could also be useful.

A. Derivation of the PDE for the average jet multiplicity
Differentiating Eq. (6.3) with respect to ξ and using (3.3), we have to NDLA But at the same level of precision we have (remembering the ξ dependence of α S ) and Differentiating with respect to E and using (3.3) again, Now to the required precision Thus we get the PDE which gives (6.5) in terms of the logarithmic variables (4.9).

B. Properties of the PDE for the average jet multiplicity
The PDE (6.17) is not so straightforward to solve explicitly. If we consider changing to variables corresponding to the sum and difference of κ and λ, we find that the PDE is separable (one obtains ODEs that are equivalent to those describing the classical mechanics of the harmonic oscillator and the quantum mechanics of the Coulomb potential), but that the boundary conditions are not. As a result, one may prove that one cannot express the boundary conditions as a Fourier-Bessel series of orthogonal functions of Sturm-Liouville type. On reflection, seeking an explicit solution in this way is perhaps ambitious, given that there is no guarantee that solutions to such a PDE can be written in terms of special functions. We thus proceed to an analysis of a rather different nature. This analysis will enable us to (i) establish that a solution exists; (ii) give a variety of infinite montonic series of upper and lower bounds on that solution; and (iii) provide explicit series solutions that enable, e.g., the asymptotic behaviour of the solution to be found in simple, closed form.
To proceed, it is useful to define x = cg 2 (2κ + µ 2 ), y = cg 2 (λ + µ 2 ), in terms of which We may recast the PDE in the integral forms where a = cgµ 4 , or 3) The first of these forms, by the way, makes clear that a physically-acceptable solution, N , if it exists, is a monotonically increasing function of x, ∀ y, and of y, ∀ x. (Proof: The mean number of jets must be ≥ 0. Thus N y > 0; N x > 0 follows by symmetry in x ↔ y.) The basic idea behind our analysis will be to first identify related PDEs for which we can find explicit solutions whose properties (such as their existence, continuity, &c.) can be checked 'by hand'. We then use these solutions as a crutch to derive properties (for example, the existence) of the solution of the original PDE.
Our first lemma is as follows. Suppose M 0 (x, y) is a solution of the PDE  ( For brevity, we denote a double integral of the type that appears on the RHS by M 0 x+y henceforth.) Proof of this follows immediately from the fact that 1 + M 0 x+y ≤ 1 + αM 0 = M 0 everywhere in D. Continuing in this vein, we define an infinite sequence whose elements are y) ∈ D and proof that this sequence is montonically decreasing follows by induction, since we have already proven that M 1 ≤ M 0 everywhere in D. We next prove that the sequence {M i } is bounded below by zero, if M 0 > 0 (which we can verify by hand given an explicit M 0 ), and hence converges. To do so, we simply note that positivity of M 0 implies positivity of M i ∀ i by induction, given (B.6) and M 0 > 0.
We now prove that the sequence {M i } converges to a solution, N of the original PDE. Indeed (B.8) where, in the penultimate step we used the Theorem of Monotone Convergence of Lebesgue integration. 2 We shall call N a supersolution of (B.1). Evidently, given any suitable α and M 0 , we can use the above results to prove that a solution exists and to find an infinite series of monotonically decreasing upper bounds on it. This raises the question of uniqueness, however: different starting points, (α, M 0 ), may lead to different supersolutions.
It is straightforward to prove the following weak version of uniqueness, for sequences that are 'nested' in the following sense. Suppose two sequences {M i } and {M i } converge on supersolutions N and N , respectively. If ∃ i, j, k, l such that elements M i ≥ M j and M k ≤ M l almost everywhere, then N = N . Proof: Unfortunately, this weak version of uniqueness does not even imply that any two supersolutions coincide, since the sequences that lead to them may be wholly contained in distinct intervals, or, even if they are nested, an element of one sequence need not exceed an element of the other almost everywhere. But this weak version is still useful where it applies, in that it may give us a more rapidly converging sequence of upper bounds on a solution.
Let us now discuss lower bounds. Similar arguments to those just given show that if we instead started with a solution, L 0 (x, y) to a PDE of the form L 0 xy = β(x, y)L 0 , (B.9) with β < 1 x+y almost everywhere on D, then we shall obtain an infinite, monotonically increasing sequence. We can, furthermore, prove that this sequence is bounded above and hence converges, provided we can show by hand that L 0 ≤ M 0 , for some M 0 that converges to a supersolution N . To wit, suppose we have shown that L 0 ≤ M 0 . It then follows by induction that L i ≤ M i , since L i+1 − M i+1 = L i −M i x+y , and thus that lim i→∞ L i ≤ N . Again, by the Theorem of Monotone Convergence, {L i } converges to a solution of (B.1) that we call a subsolution and denote by N . By a straightforward generalization of the result just proven, a subsolution is less than or equal to any supersolution for which it can be shown that any one element in the sequence defining the former is less than or equal to (almost everywhere) any element in the sequence defining the latter.
Yet again, this result is insufficient to establish uniqueness, but it does provide a way to obtain both upper and lower bounds on a given solution. Obviously, if we can prove (or assume) uniqueness, then all of the aforementioned super-and subsolutions coincide, meaning that any of the aforementioned sequences can be used to bound the solution to arbitrary, known precision.
We can prove uniqueness of the solution heuristically out to any finite (x, y) by explicit construction. Given N = 1 on the boundary, we divide the intervals [a, x] and [a, y] into regions of size δx and δy and find the unique solution N (a + δx, y) 1+δx y a dy a+y , with a similar result for N (x, a+δy). We then extrapolate to the region [a + δx, a + 2δx] and so on. Given that we have proven the existence of the solution, we do not need to worry that the extrapolations in the x and y directions might not coincide. To make this proof rigorous outside a neighbourhood of the boundary, we would then have to take the limit δx, δy → 0. We will content ourselves with assuming that the solution is unique.
As we stressed above, all of our results on super-and sub-solutions are contingent on showing the existence of solutions to the equations (B.4) and (B.9), for suitable α and β, satisfying the BCs, and showing that they have the desired properties. We do this by supplying explicit solutions for functions α, β of the combination (x−a)(y−a), for which the PDE reduces to an ODE, and for which the boundary conditions reduce to a single boundary condition at (x − a)(y − a) = 0.
To obtain lower bounds, once may start, for example, with the solution L 0 = 1 of the PDE L 0 xy = 0, since 1 x+y > 0 everywhere in D. We thus obtain L 1 ≥ 1 + (x + y) ln(x + y) − (x + a) ln(x + a) − (a + y) ln(a + y) + 2a ln 2a , (B.10) with subsequent L i given in terms of polylogarithms. At fixed y and large x, for example, we deduce that L 1 ∼ (y − a) ln x. More generally, L n ∼ ((y−a) ln x) n For a better lower bound, we can, provided a > 1, solve L 0 xy = 2 (x + y) 2 L 0 , (B.13) whose solution is L 0 = xy + a 2 a(x + y) . (B.14) It is soothing to verify explicitly that the resulting subsolution has the same asymptotic behaviour, at large x and fixed y, as the subsolution starting from L 0 = 1 derived above. Indeed, in this limit, L 0 ∼ y a , L 1 ∼ y 2 −a 2 2a log x, and so on, with L n ∼ (y−a) n (y+na) a(n+1)! log x n! . At large n, we thus find L n ∼ (y−a) n n! log x n! , exactly as before.