The unified transform for evolution equations on the half‐line with time‐periodic boundary conditions *

This paper elaborates on a new approach for solving the generalized Dirichlet‐to‐Neumann map, in the large time limit, for linear evolution PDEs formulated on the half‐line with time‐periodic boundary conditions. First, by employing the unified transform (also known as the Fokas method) it can be shown that the solution becomes time‐periodic for large t . Second, it is shown that the coefficients of the Fourier series of the unknown boundary values can be determined explicitly in terms of the coefficients of the Fourier series of the given boundary data in a very simple, algebraic way. This approach is illustrated for second‐order linear evolution equations and also for linear evolution equations containing spatial derivatives of arbitrary order. The simple and explicit determination of the unknown boundary values is based on the “ Q ‐equation”, which for the linearized nonlinear Schrödinger equation is the linear limit of the quadratic Q ‐equation introduced by Lenells and Fokas [Proc. R. Soc. A, 471, 2015]. Regarding the latter equation, it is also shown here that it provides a very simple, algebraic way for rederiving the remarkable results of Boutet de Monvel, Kotlyarov, and Shepelsky [Int. Math. Res. Not. issue 3, 2009] for the particular boundary condition of a single exponential.


INTRODUCTION
The unified transform, also known as the Fokas method, was introduced by the first author in 1997. 1 It provides a novel approach for solving boundary value problems for linear and integrable nonlinear partial differential equations (PDEs). For integrable PDEs the unified transform is the proper generalization from initial to boundary value problems of the famous inverse scattering transform (IST). Taking into consideration the seminal contributions of Harvey Segur in the development and applications of the IST, we hope that our paper fits well in this issue.
The unified transform has been employed by many authors for the investigation of evolution PDEs in one space variable. For such equations, among the most important results obtained via this method are the following: (i) Linear equations formulated on the half-line or a finite interval have been analyzed by Deconinck, Fokas, Pelloni, and collaborators. [2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17] (ii) Numerical techniques for linear equations are developed in Refs. 18-22. (iii) Novel results in spectral theory are derived in Refs. 23-25. (iv) Evolution PDEs defined on moving domains are analyzed in Refs. 26-28. (v) Unexpected results in the area of null controllability are presented in Ref. 29 (vi) Integrable nonlinear evolution equations have been analyzed by many authors, see, for example, Refs. 30-38. (vii) Nonlinear integrable equations with the so-called linearizable boundary conditions are investigated in Refs. [39][40][41][42][43]. For this type of special boundary conditions, the unified transform yields a solution, which is as effective as the solution of the corresponding initial value problem on the line obtained via the IST. (viii) It is shown in Refs. 44, 45 that -periodic initial conditions belong to the linearizable class. Hence, remarkably, such problems can be solved as effectively as the usual initial value problem on the full line. (ix) Himonas, Mantzavinos, and Fokas have initiated a program of study for using the unified transform to derive results regarding the well-posedness of arbitrary nonlinear evolution PDEs. [46][47][48][49][50][51] It is well known that for the analysis of elliptic PDEs, one has to face the difficulty of analyzing the Dirichlet-to-Neumann map. Actually, this problem also arises in the analysis of boundary value problems for evolution PDEs. Consider, for example, the case of an evolution equation involving spatial derivatives of order up to two, formulated on the half-line. Let 0 and 1 denote the Dirichlet and Neumann boundary values at = 0, i.e., For a well-posed problem, either 0 or 1 or a relationship between them is prescribed as a boundary condition. This corresponds to a Dirichlet or Neumann or Robin boundary value problem, respectively. The generalized Dirichlet-to-Neumann map involves the determination of the unknown boundary value in terms of the given initial and boundary conditions. A vital ingredient of the Fokas method is the so-called global relation. For linear evolution equations, this is a linear algebraic equation coupling appropriate integral transforms of all boundary values. By utilizing the global relation it is possible to determine the generalized Dirichlet-to-Neumann map for a linear evolution equation whose highest spatial derivative is of arbitrary order. However, for nonlinear integrable evolution equations the global relation is nonlinear. Thus, the generalized Dirichlet-to-Neumann map is characterized, in general, in terms of a nonlinear formalism. Despite this serious difficulty, substantial progress has been made via the usage of the global relation: First, it was realized in Ref. 31 (2), and the corresponding Neumann boundary value also asymptotes to a single exponential: (iii) Rigorous well-posedness results for the Dirichlet-to-Neumann map for decaying Dirichlet data are derived in Ref. 56. (iv) For vanishing initial data, if the Dirichlet boundary condition is a sine wave, i.e., (0, ) = sin , ∈ ℝ, then the Neumann boundary function (0, ) for the NLS, 34 and the Neumann boundary functions (0, ) and (0, ) for the modified Kortewegde Vries equation, 57 can be computed up to and including terms of ( 3 ) and, at least up to this order, are asymptotically time-periodic for → ∞. A similar result has been obtained for the KdV equation, 58 where it is shown that when the Dirichlet condition is a sine wave, then the Neumann boundary values are asymptotically time-periodic at least up to second order in perturbation theory. This perturbative approach introduced in Ref. 34 involves heavy calculations, making it virtually impossible to go further than terms of ( 3 ), but gives however a strong indication that (asymptotically) time-periodic Dirichlet boundary conditions lead to asymptotically time-periodic Neumann boundary values.
It is known that, for integrable evolution equations formulated on the full line, the power of the integrability formalism becomes evident in the asymptotic analysis for large . This is also true for boundary value problems analyzed via the Fokas method. Indeed, for boundary conditions, which vanish for large , the unified transform gives rise to a Riemann-Hilbert formulation with explicit and dependence. Hence, it is possible to determine the structure of the solution without determining the Dirichlet-to-Neumann map. 59 As a result of the fact that this map remains undetermined, certain constants in the relevant asymptotic formulas remain unknown.
Regarding the case of -periodic boundary conditions, a new approach for the large asymptotic analysis was introduced in Ref. 60 (see also Ref. 58). This approach is based on the investigation of a certain quadratic equation. We will refer to this equation, which can be obtained directly from the Lax pair evaluated at = 0, as the -equation. It was shown in Ref. 60 that for a Dirichlet boundary condition, which asymptotes to a -periodic function for large that can be expanded in a power series of , it is possible to obtain the Neumann boundary value via a series expansion in , which is explicitly determined in terms of the Dirichlet data.
In the present paper, it is shown that the -approach is truly powerful: (i) For the case of the linear version of the NLS, namely, for the equation it was shown in Ref. 60 that, if the Dirichlet condition is asymptotically -periodic, then for large the Neumann boundary value is also -periodic. Moreover, an explicit formula was obtained in Ref. 60 relating the coefficients of the two Fourier series characterizing the asymptotic Dirichlet and Neumann values. It is shown in Section 2 that the linear version of the -formalism provides a much simpler way for computing asymptotically the Dirichlet-to-Neumann map of Equation (4). Similar results are obtained in Section 2 for the heat equation, and for the diffusion-convection equation, (ii) For large , this new approach provides, in a straightforward manner, the asymptotic form of the generalized Dirichlet-to-Neumann map of linear evolution equations containingderivatives of arbitrary order, with asymptotically time-periodic boundary data and with a sufficiently decaying initial condition. We illustrate the general approach in Section 3 using the example of the Stokes equation, The specific implementation of the new method to Equations (4)-(7) makes clear how this method can be easily applied to any linear evolution equation with constant coefficients. (iii) The determination of the large -asymptotics of the Dirichlet-to-Neumann map of the diffusion-convection equation (6), which contains the heat equation as a special case for = 0, is presented in Section 4 via the unified transform method. The derivation follows similar steps to those used in Ref. 60. However, it presents a slight generalization of the analogous result of Ref. 60, because here we derive the large asymptotics of ( , ) as opposed to (0, ). Using almost identical steps it is possible to derive an expression for the large behavior of ( , ), which shows that ( , ) becomes -periodic. By evaluating the expression for ( , ) derived in this way at = 0, we obtain explicitly the relationship between the Fourier coefficients of the Dirichlet and the Neumann values. By comparing this derivation with the very simple procedure of the -approach, the advantage of the latter becomes evident. Analogous computations for the Stokes equation (7)

THE LARGE BEHAVIOR OF THE DIRICHLET-TO-NEUMANN MAP FOR SECOND-ORDER LINEAR EVOLUTION EQUATIONS WITH -PERIODIC BOUNDARY CONDITIONS
The derivation of the -equation is based on the Lax pair formulation. In turn, this is related to the first step of the Fokas method, namely, rewriting a given linear PDE in a divergence form. This can be achieved in a variety of ways, including the use of the adjoint. Let ( , ) satisfy the second-order evolution PDE where 0 , 1 , and 2 are complex constants. We assume that ( , ) is sufficiently smooth and that both ( , ) and its -derivatives decay sufficiently fast as → ∞ for each ≥ 0.
The formal adjoint can be obtained from Equation (8) by replacing ∕ and ∕ with − ∕ and − ∕ , respectively. Hence, Multiplying Equations (8) and (9) by and , respectively, and then adding the resulting equations, we find Equation (8) We assume that ReΩ( ) ≥ 0 for real, so that the initial value problem is well-posed. 62 We observe that, e − +Ω( ) is a solution of Equation (9). Replacing, in Equation (10), by this exponential, we find the following one-parameter family of divergence forms: This equation implies the existence of the function e − +Ω( ) ( , , ), where Simplifying, we obtain the following Lax pair of Equation (8):

The -formulation
Suppose that Equation (8) is valid in a given domain . Equation (12) and Green's Theorem imply that where denotes the boundary of . In the particular case that Equation (8) is formulated on the half-line, Equation (15) becomes the so-called Global Relation: Let ( , ) satisfy the -dependent part of the associated Lax pair evaluated at = 0, i.e., ( , ) + Ω( ) ( , ) = 2 1 ( ) + ( 2 + 1 ) 0 ( ), > 0, ∈ ℂ. Then, Comparing this equation with the Global Relation (17), it follows that if we choose For the integrals appearing on the right-hand side of Equation (19), we have − ≥ 0. Let us assume that 0 and 1 are bounded functions. Thus, as → ∞, the second and third terms of the right-hand side of Equation (19) are bounded and analytic in the domain Moreover, ( , 0) is bounded and analytic for Im ≤ 0. Thus, ( , ) is well defined in The unified transform yields a representation, which involves the curve − , which is the boundary of the domain − defined by We observe that the domain − is the complement (in the lower half of the complex -plane) of the domain (23) of validity of . However, Equation (21) implies that, actually, ( , ) has a much larger domain of analyticity. Namely, this equation extends the domain of validity to Im ≤ 0 (the entire lower half of the complex -plane including the real axis), which encompasses − , i.e., − and its boundary. Thus, we will require that ( , ) is free of singularities in − . Let us assume that for large both 0 and 1 asymptote sufficiently fast toward smooth timeperiodic functions of the general form where (0) , (1) ∈ ℂ, ∈ ℤ and > 0. We seek a solution of Equation (18) of the form Substituting the above asymptotic expressions for 0 , 1 , and in Equation (18), we find Remarkably, the condition that ( ) does not have poles for ∈ − , clearly imposes a relationship between (1) and (0) . Actually, if the denominator of Equation (27) becomes zero for some 's then, because is purely imaginary, those 's should lie on the contour defined by ReΩ( ) = 0. Hence, any singularities of ( ) in − occur in fact on its boundary − . In other words, the condition that the singularities of ( ) on − are removable, determines the Dirichlet-to-Neumann map.
Remark 1. The assumption (25), i.e., that (asymptotically) -periodic Dirichlet boundary conditions lead to asymptotically -periodic Neumann values, can be shown to hold for appropriate initial-boundary value problems using the Fokas method and the method of steepest descent, see Section 4 for details. In particular, to avoid problems with the steepest descent arguments, we have to make sure that the denominator involved, which for linear evolution PDEs is of the form Ω( ) + , does not become zero on the steepest descent contour. The corresponding steepest descent contour is the contour that passes through a saddle point of Ω( ) and on which the imaginary part of Ω( ) is constant. So, we can a priori consider appropriate -periodic boundary conditions to avoid problems with the steepest descent arguments.
In what follows we consider three particular examples.

2.2
The linearized NLS equation Let ( , ) solve the linearized version of the NLS, namely, In this case, 2 = , 1 = 0 = 0, and we assume that Denoting = + with , ∈ ℝ, the domain − defined by Equation (24) takes the following form: which is the third quadrant of the complex -plane. The singularities of ( ) occurring on − are given by The condition that the above singularities are removable, implies the following relationship between (0) and (1) : with ( ) defined by Equation (31). This relationship was first given in Ref. 60.

The heat equation
Let ( , ) solve the heat equation, i.e., In this case, 2 = 1, 1 = 0 = 0 and we assume that The domain − , defined by Equation (24), takes the following form (describing a wedge-shaped domain in the lower half of the complex -plane): where = + with , ∈ ℝ. The condition that the singularities ( ) of ( ) occurring on the boundary are removable, implies the following relationship between (0) and (1) : with ( ) defined by Equation (36).

The diffusion-convection equation
Let ( , ) solve the diffusion-convection equation: In this case, 2 = 1, 1 = , and 0 = 0 and we assume that The singularities of the denominator lie at For the diffusion-convection equation, the domain − takes the form where = + with , ∈ ℝ. It can be readily checked that the singularities that lie on − are given by the solutions with the minus sign of Equation (40), which we denote by 2 ( ). Hence, the condition that these singularities are removable, implies the following relation between and (1) : with 2 ( ) defined by the solutions with the minus sign given in Equation (40).

THE LARGE BEHAVIOR OF THE DIRICHLET-TO-NEUMANN MAP FOR AN ARBITRARY LINEAR EVOLUTION EQUATION WITH -PERIODIC BOUNDARY CONDITIONS
Let ( , ) satisfy the following general evolution PDE on the half-line, which has spatial derivatives of arbitrary order: where Ω( ) is a polynomial of the complex variable of degree and ReΩ( ) ≥ 0 for real. We assume that ( , ) is sufficiently smooth, up to and including the boundary, and that both ( , ) and its -derivatives decay sufficiently fast as → ∞ for each ≥ 0.
A particular solution of Equation (43) is It is shown in Ref. 62 that Equation (43) can be written in the form where { ( )} −1 0 can be calculated explicitly in terms of Ω( ) as follows: Hence, the associated Lax pair is As before, letting ( , ) satisfy the -dependent part of the associated Lax pair evaluated at = 0, we find where the 's are defined by ( ) ∶= (0, ). Assuming the following expansions for large , and substituting the above expansions, along with Equation (26), into Equation (48), we find In the same way as in Section 2, we can show that the singularities of ( ) on − have to be removable. Hence, this condition determines the Dirichlet-to-Neumann correspondence also for the general linear evolution PDE (43). We illustrate the general approach using the example of the Stokes equation.
As before, we require the singularities of ( ) on − to be removable, where for the Stokes We observe that − = 2 ∪ 3 , see Figure 1. Let us denote with 1 ( ), 2 ( ), and 3 ( ) the zeros of the denominator, which lie on the boundaries 1 , 2 , and 3 , respectively. Thus, the singularities 2 ( ) and 3 ( ) are the ones that lie on − . Because these two singularities should be removable, we need the numerator of (53) to vanish for = 2 ( ) and = 3 ( ), so we arrive at the following system for the coefficients (0) , (1) , and (2) : Let us assume that we consider the Dirichlet boundary problem of the Stokes equation and, hence, the coefficients (0) , ∈ ℤ, are given. We can then proceed to solve the above system for the unknown Neumann coefficients (1) and (2) in terms of the Dirichlet coefficients (0) . From Substituting this expression into Equation (56b), we find the first set of Neumann coefficients (1) (in terms of the Dirichlet coefficients): where we have used the Vieta formula 1 ( ) + 2 ( ) + 3 ( ) = 0. From Equation (56b) we have Substituting the value of (1) from (58) into Equations (57) and (59) and adding the resulting equations, we find where we have again used that 1 ( ) + 2 ( ) + 3 ( ) = 0. Hence, the second set of Neumann coefficients (2) (in terms of the Dirichlet coefficients) is given by Thus, we have determined the Dirichlet-to-Neumann map for the Stokes equation.

THE UNIFIED TRANSFORM FOR THE LARGE ASYMPTOTICS OF THE DIFFUSION-CONVECTION EQUATION Proposition 1. Let ( , ) satisfy the following Dirichlet problem for the diffusion-convection equation (6) on the half-line:
where the initial and boundary data are compatible at the origin. The solution ( , ) of this problem is assumed to be sufficiently smooth (up to and including the boundary) and both ( , ) and its -derivatives are assumed to decay sufficiently fast to zero as → ∞ for each ≥ 0. The initial condition 0 ( ) is therefore taken to have sufficient smoothness, and both 0 ( ) and its derivatives are taken to decay sufficiently fast to zero as → ∞. We further suppose that the Dirichlet boundary condition 0 ( ) is bounded, sufficiently smooth, and for large asymptotes sufficiently fast toward a periodic function where > 0, the prime denotes that with (1) where the 2 ( )'s, ∈ ℤ, denote the roots with the minus sign given in Equation (40).
where the contour + is depicted in Figure 2 Because we are interested in finding the Neumann boundary value in terms of the Dirichlet boundary condition, we differentiate both sides of Equation (70) with respect to and we find In what follows, we will analyze separately the terms of Equation (71). Let us begin by analyzing the first two integrals appearing in Equation (71). Replacing with − + in the second of these integrals, and using the definition of̂0 from Equation (68), the sum of these integrals takes the form In the third step of the above computation we integrated by parts, so that in the next step we can now change the order of integration and obtain well-defined -integrals. Changing the order of integration in Equation (72) and completing the squares in the exponents, we find that Equation (72) can be written as (75) The above expression vanishes as → ∞, by also taking into account that 0 ( ) anḋ0( ) are assumed to have sufficient decay as → ∞, see the assumptions below (62). Thus, for → ∞ the -derivative of , given by Equation (71), reduces to where where we can differentiate under the integral sign by Leibniz's rule, because we assume that 0 is bounded and we consider strictly positive . In what follows, to avoid convergence issues, we assume that > 0 and only at the end of this section we will take the limit → 0 + to arrive at the Neumann boundary value. Assuming that in the large time limit 0 ( ) tends to a periodic function sufficiently fast, where > 0 and the prime denotes that Let us rigorously justify the above formula. From Equation (78) we have that for any > 0, there exists a ( ) > 0, such that ≤ , for all > ( ).
Note that Equation (79) can be justified in a similar way for = 0, i.e., the heat equation, but in this case the convenient contour to which we deform + is taken to be where = + with , ∈ ℝ. For completeness, we note that for other equations (such as the linearized NLS equation) more caution may be needed for the justification of Equation (79).
From Equation (79) we find where in the first line we were able to interchange the order of integration and summation because of the assumption in Equation (64).
The singularities of the expression in Equation (89) lie at = 1,2 ( ), ∈ ℤ * , as defined in Equation (40). We observe that the roots 1 ( ) with the plus sign of (40), i.e., lie on + for all ≠ 0. However, they are removable singularities, because the following limit exists: Thus, we can deform the contour + tõ+, which passes below the singularities = 1 ( ), ∈ ℤ * , see Figure 2(B). We split the integral appearing in Equation (89) into two terms, which we call 1 ( , ) and 2 ( , ), as follows: Note that we were able to interchange the order of integration and summation in the term 1 ( , ) of Equation (92) because of the absolute convergence assumption on the Fourier coefficients (0) , see Equation (64). We first focus on 2 ( , ) and show that it vanishes as → ∞. The term 2 ( , ) involves an integral of the form where Its large -asymptotics (for each given ) can be computed via the method of steepest descent. 63 We deform the contour̃+ to a new contour ′ on which Φ( ) has a constant imaginary part and which passes through the saddle point 0 = 2 of Φ( ). Thus, the contour ′ is the horizontal line going through 0 = 2 . Observing that 0 = 2 is a simple saddle point of Φ, meaning that Φ ′ ( 0 ) = 0 and Φ ′′ ( 0 ) ≠ 0, we find , as → 2 . (95) It can easily be checked that ( . Without loss of generality, the latter expression can be assumed (for > 0) to converge to a nonzero value. Even if the (0) 's would happen to be such that ( 2 , ) converges to zero, we could proceed to compute higher order -derivatives of ( , ) at 0 = 2 until one of them would be nonzero; in this case we would get an even faster rate of decay in Equation (97) as → ∞. Hence, assuming that ( 2 , ) ≠ 0, we have Using the formula given by eq. (6.4.9) of Ref. 63, we find Note that in the above analysis we have used the absolute convergence assumption on the Fourier coefficients (0) , Equation (64), to assure convergence and to find the derivative of the infinite sum. Also, note that for the special case = 0, which corresponds to the heat equation, we find = 2, = 3, 0 = 0, and 0 ( ) = 1 2 2 (0, ). Substituting these values in Equation (97), and noting that Φ(0) = 0, we again find that 2 ( , ) vanishes as → ∞. So, from (92) and (98) we have Hence, recalling Equation (76), we observe that ( , ) is asymptotically -periodic. In an almost identical way, one can find an expression for the large behavior of Equation (70) and show that ( , ) itself also becomes -periodic for → ∞.
At this point, it is important to observe that ( ) is bounded at = Γ 2 , where > 0 and Γ ≥ 0. This can be seen by substituting the asymptotic forms of 0 ( ) and 1 ( ) from Equation (116) into the -equation (115), because one then finds ( Γ 2 , ) = exp[ (2 2 − Γ 2 ∕ 2 ) ], ∈ ℝ, when → ∞. So we have to require that 4 2 − 2 2 + = 0 when the denominator is zero, i.e., when = Γ 2 . This can also be inferred from Equation (121), where for = Γ and = Γ 2 (with Γ real and nonnegative) only the second term survives, which then necessarily must be equal to zero. Under the assumption that ( . Hence Γ = √ 2 2 − , and thus we find in agreement with Equation (3b). It may be noted that the inequality ≤ −6 2 of Boutet de Monvel, Kotlyarov, and Shepelsky is not reproduced, yet the present analysis provides an elementary and straightforward approach for determining the large behavior of the Dirichletto-Neumann map in the case of periodic single-exponential boundary functions.
Boutet de Monvel et al. 55 showed that the coefficient of the asymptotically periodic Neumann boundary value is as presented in Equation (3), i.e., is either real or purely imaginary if is outside the range (−6 2 , 2 ). Inside the range ∈ (−6 2 , 2 ) the asymptotic form of the Neumann boundary value is not of the simple form e , 64 so then Equation (116) would not be a valid substitution to make. Hence, the above analysis concerning real and purely imaginary (but not more general complex values of ) covers all cases of interest for the focusing NLS with boundary data that asymptote to periodic single exponentials.

CONCLUSIONS
The goal of this work is to elucidate the effectiveness of the use of the -formulation. For integrable nonlinear evolution equations, the -equation was introduced in Ref. 60 in connection with the NLS. In the present paper we have established the following: (i) For linear evolution equations, the analysis of the -equation yields the large asymptotics form of the generalized Dirichletto-Neumann map for -periodic boundary conditions, in a very simple, algebraic way. (ii) For the NLS, it reproduces the remarkable results of Ref. 55, again in a simple, algebraic manner. We expect our method to give useful results also for other integrable nonlinear PDEs, which are amenable to the Fokas method. For linear evolution equations, the -equation is the -part of the Lax pair evaluated at = 0. It is interesting to note that when the unified transform was first introduced, it was implemented via the Lax pair formulation. Later, it was realized that (in the linear case) it could be derived in a straightforward manner, avoiding the Lax pair connection. Of course, the Lax pair approach remains indispensable for the application of the unified transform to nonlinear integrable PDEs. The results presented in this paper clearly show that, even for linear evolution equations, the connection with the Lax pair is very useful. Indeed, in terms of simplicity, the -approach (which is based entirely on the Lax pair formulation) provides the most unexpected results, so far, of the Fokas method. The effectiveness of this new approach becomes apparent by comparing it with the lengthy derivation in Section 4 and earlier implementations of the unified transform in Refs. 55, 60. Taking into consideration that the classical investigations of linear PDEs did not take into account that linear PDEs admit a Lax pair formulation, 65 it is not surprising that the -approach was missed in these previous investigations.
The -approach is based on the assumption that the unknown boundary values become periodic for large . Although it is rather complicated to determine their precise asymptotic form, it is often much simpler to show that they indeed become -periodic as tends to infinity. Using the general formulation of the Fokas method and the standard asymptotic technique of steepest descent, it can be shown in many instances that, if the given boundary conditions areperiodic, the unknown boundary values also become -periodic as tends to infinity. Then, the -approach provides a most effective way for computing explicitly the Fourier series coefficients of the unknown periodic functions in terms of the Fourier series coefficients of the given boundary conditions.

A C K N O W L E D G M E N T S
A. S. Fokas is supported by the EPSRC in the form of a senior fellowship. M. C. van der Weele gratefully acknowledges support from the Onassis Foundation (Scholarship ID: F ZQ 004-1/2020-2021) and from the Foundation for Education and European Culture, founded by Nicos and Lydia Tricha. During the early stages of this work she was also supported by the Cambridge Trust and Christ's College via the "Vice-Chancellor's and Christ's College Warwick Scholarship," as well as by The A. G. Leventis Foundation, for which she is very thankful. Both authors want to thank the anonymous referees for very insightful comments and suggestions.

D ATA AVA I L A B I L I T Y S TAT E M E N T
Data sharing not applicable to this article as no data sets were generated or analyzed during the current study. This article describes entirely mathematical/theoretical research and does not rely on any experimental data.