One of the important quantities one can calculate using the mean-field method for the polynomial models is the position of the critical curve within the parameter plane of each model. This is done in section 3 of chapter 2 of the book [1], for the case of the model. The result can be extended to the models, and in this general case the mean-field result for the critical curve of the polynomial model is:

which determines implicitly. Since it is not possible to solve this equation analytically in order to obtain in explicit form, we will be concerned here with the numerical solution of this equation, so that we may plot the graph of this function. The integrals can be written in terms of the parabolic cylinder functions [3],

However, this does not help us to solve the equation numerically because these functions are not readily available in a computer implementation. Note that for even the integrals can be written as integrals on a variable ,

Hence it is not surprising that for even they can be written in terms of the error function, which can in fact be done through the functions, with integer . The functions can be written in terms of the error function , which is given by [4]

For the first few functions we have [5]

From the first two one can get all the functions in this sequence by repeated use of the recurrence relation [6]

which can also be written in the form

Using all this one can solve the equation numerically for the case of
even , one case at a time, using the error functions `erf` or
`derf`, which are readily available for use with the `g77`
Fortran compiler. For example, for
, which is the simplest case,
the equation is

Given a value of , one can solve this equation numerically for , and vice-versa. However, this solves only the even- cases, not the odd- cases, including .

For all values of one can make a change of variables in order to simplify the original integrals in a way that is convenient for numerical purposes. Since the solution for is obviously , we can assume that and consider the variable given by

in terms of which the equation becomes

Defining now a new parameter given by

we may write the equation in the form

or, multiplying above and below by ,

Given a value of , this allows one to calculate using a single pair of numerical integrations, and afterward to calculate the corresponding from , and . The parameter is a new parameter for the critical curve. In order to determine its range of variation as we travel along the curve from end to end, consider first the asymptotic line of the critical curve, for large values of and , which can be obtained through the representation of the equation in terms of the parabolic cylinder functions [7], and is given by

where is the critical point of the sigma model over the arc at infinity. In this case we have for the expression

Therefore corresponds to and hence to the sigma-model critical point over the arc at infinity, and when

If we now consider the tangent to the critical curve at the Gaussian point, which can be obtained by differentiating equation (1) implicitly with respect to and , then expressing the resulting integrals in terms of the function [8], and which is given by

we see that in this case we get for the expression

Therefore corresponds to and hence to the Gaussian point, and when

In short, we have the following limiting behaviors:

We have therefore the following general algorithm for the numerical solution of the equation:

- Choose a value for in
.
- Calculate
by numerical integration,

- Calculate the corresponding by the formula

The resulting pair is a point of the critical curve of the model, in dimensions. Probably the best way to choose values of so that the resulting points are distributed fairly homogeneously along the critical curve is to use the equation of the tangent line at the Gaussian point, choosing some , integer values of , and using

Since the functions being integrated are centered around a single maximum that moves to large values of the variable when is large, in order to integrate them efficiently we must know the location of the points of maximum of the functions. The basic idea for the integration of one of these functions is to start at the point of maximum of the function, and then to integrate to both sides, down towards zero and up towards infinity. Let us consider then the functions

where is an integer. Except for the case this is zero for , and in all cases for . Taking the derivative we get

Besides vanishing at (for ) and at , which are clearly points of minimum, this can be zero for

and since only the sign can be kept, and we get

This is correct in all cases including . For we have that, if , then , and if , then . Note that the only case in which the maximum is at zero is with , in all other cases the maximum is at a positive and non-vanishing value of .

The functions examined above have the form of a pulse around their points of maximum, but this pulse has a variable width, which decreases as the parameter increases to positive values. This makes is harder to find an appropriate numerical integration interval for them. For the cases we can perform one more transformation of the integrals, in order to improve this situation. We make a transformation of integration variables to the variable and hence write

where the integrand has now a fairly constant width around its maximum. In this case the location of the point of maximum is a bit different, that is, is not exactly . The integrand is now

and its derivative is given by

Note that, due to the possibly negative powers of , the cases and should be examined separately. For we have

which is zero at the point of minimum at infinity and at , which is therefore the point of maximum. Since this is valid only for , for the function is monotonically decreasing with its maximum at . For we have

Note that in this case the function has an infinite derivative at the point , which makes it unsuitable for numerical integration. Its point of maximum can be obtained as one of the roots of a quadratic polynomial and is in fact included in the general case, which we proceed to analyze. If we multiply it by , the derivative of the general case becomes

Except for the cases and this vanishes at the points of minimum at and at , and for all cases it vanishes at the point of maximum given by

where since only the sign can be kept, and we get

This is in fact correct in all cases including and . For we see that this gives the correct answer: if , then , and if , then . Note that the only case in which the maximum is at zero is with , in all other cases the maximum is at a positive and non-vanishing value of .