Molecular response for nematic superconducting media in a hollow cylinder: a numerical approach

In the context of the Ginzburg–Landau formalism proposed by Barci et al. in 2016 for nematic superconductivity, and by performing a numerical treatment based on the shooting method, we analyze the behaviour of the radial distribution of the nematic order parameter when the superconducting order parameter reaches the typical non-trivial saddle point. We consider the cases of a hollow cylindrical domain, with a disk or an annular domain as its cross section, whether the order parameter is subjected to Neumann or Dirichlet boundary conditions. We conclude that depending on the corresponding situation a non trivial solution holds if certain relations between the radii are satisfied. Moreover, we observe a saturation effect on each instance that constitutes a purely geometrical consequence on the relation between the typical sizes and shapes of the samples. These conclusions can be useful for further experimental realizations and extensions to the interaction of the nematic (superconducting) order parameters with electromagnetic fields.


Physical Preliminaries
Critical phenomena have been extensively studied to determine the behavior of several properties as a function of the temperature.From a macroscopic point of view and in the low temperature limit, Ginzburg Landau equations have been successfully applied to explain superfluidity, superconductivity [8,11,13,18,21] and liquid crystals [7,16], among others.In the context of superconductivity, the decrease of a magnetic field due to the formation of Cooper pairs, the mathematical formalism of the Ginzburg Landau theory have also been amazingly developed [17] to understand the properties of the energy functional and its minimizers, and as a complement to explain the interaction between superconducting vortices.The latter stems from the original problem proposed by Abrikosov for unbounded domains in 1957 [1].In this sense, the attractive (repulsive) nature of two equal type I (II) superconducting vortices was explained theoretically [12] and experimentally [19].Besides, the repulsive interaction between two equal type II vortices was demonstrated in bounded domains by applying an electromagnetic analogue of two mutually inducting coils in a long hollow cylinder domain [15].
On the other hand, liquid crystals and particularly nematic ones, where the molecules are oriented with respect to a given vector, have gained attention due to their exotic behaviour in presence of an electromagnetic field.In this context, the analysis of optical effects such as possible Fréedericksz transitions [9] and the formation of solitons [4,5] could influence the overall response of nematic liquid crystals in any of their multiple biological and electrical applications.
Starting from the simplest PDW state, Barci et al. proposed in 2016 [3] a Ginzburg Landau model for nematic superconductivity that is characterized by two order parameters, which shows the translational and rotational symmetry breaking of the system.This model could represent a simple platform to elucidate the optical response of these and similar materials and to compare with experiments [20].In fact, there is already theoretical works that support the modification of the usual Fréedericksz threshold in purely nematic liquid crystals, when the system is considered as a nematic inhomogeneous superfluid in presence of an electric field [10].
In this article, we perform a numerical analysis of an approximation for the dimensionless 2D Ginzburg Landau equation for the nematic order parameter of a nematic superconductor media in a hollow cylinder without electromagnetic fields, which includes both geometries: a disk and a ring.Due to translational symmetry [3], the nematic order parameter has the form G( x) = T ( x) exp 2iα( x), which T 2 representing the mean value of nematic molecules and α is its phase.In addition, the superconducting order parameter ψ = F 0 exp iθ is such that F 0 = 1 and the phase θ is quantized in order to obtain A − ∇θ = 0, B = 0; here, A is the magnetic vector potential and B the magnetic field such that B = ∇ × A.
In this framework, it is natural to consider Dirichlet (D) boundary conditions (b.c.); however, we also consider Neumann (N) b.c., due to its relevance in the context of superconductor models.Therefore, we pose the problem on a ring R i x R e and we study the behaviour of T , the radial distribution of the order parameter, in a wide range of regimes (including both (D) or (N) b.c.): from the (outer) threshold R e ∼ R + c to R e 1 and, for fixed R e , from R i ∼ 0 to the (inner) threshold R i ∼ R a (R e ) − ; the existence of these thresholds being one of the main goals of this paper, see Theorem 1.It is also worth to mention that the radial variable r is scaled with the coherence length of the system and so are the inner and outer radii R i and R e .
Regarding the numerical treatment, a broad type of results are available in the literature, most of them coming from superconducting models, where the parameter is given by the penetration length which is assumed to be small, see [2,14] and references therein.However, our treatment relies on an easy-to-use algorithm based on the shooting method performed with high-precision arithmetic that allows us to easily access the qualitative analysis without any further restriction on the size of the samples.In particular we explore and describe the behaviour of the solutions for values near the reported thresholds; besides, we also explore the saturation effect present near the boundaries of the domain.
This paper is organized as follows: first, we formulate the problem and present some numerical issues mainly concerned with the saturation effect exhibited by the solutions, that is the existence of small intervals in which some kind of transition takes place; then, we restrict to the problem posed in a disk with (D) b.c. and describe in full detail the proposed algorithm; next, we use this algorithm to characterize the solutions (for the problem posed in a disk with (D) b.c.) under qualitative different regimes focusing mainly in the saturation effect exhibited by the solutions and we provide estimates for the length of these transitions as a function of the radius; thus, we take profit of the results presented and make the (minor) changes to adapt the algorithm for the remaining situations: we first consider the disk with (N) b.c. and then we move to the ring with (D) and (N) b.c., in this order.Thus, we make a rigorous approach to establish the existence of thresholds.Finally, we discuss the results and its physical relevance.

Formulation of the problem and notation
We are interested in the behaviour for the solutions of the Ginzburg Landau equation Without loss of generality and in systems where a rotational symmetry is broken, a suitable solution for Eq.( 1) has the form G( x) = T ( x) exp 2iϕ, where ϕ ∈ [0, 2π] is the rotational angle.Therefore, we see that previous equation is restricted to the invariant subspace S 2 = span{J 2 (λ k r) exp i2ϕ : k = 1, 2, . ..}; here (r, ϕ) are polar coordinates, J 2 is the regular Bessel function of order two and the eigenvalues satisfies J 2 (λ k ) = 0 or J 2 (λ k ) = 0, depending whether Dirichlet (D) or Neumann (N) boundary conditions (b.c.) are used.Introducing polar coordinates, we write previous equation as follows: with 0 R i < R e .Here, the variable t is normalized by the coherence length of the system, which help us to measure the saturation from T = 0 to T = 1.Since T is related to the square root of the mean density of the nematic molecules, it is natural to consider Dirichlet boundary conditions.However, regarding the superconducting framework in which the natural boundary condition is given on the radial derivative, and since from a mathematical point of view an unified treatment is possible, both Dirichlet and Neumann boundary conditions will be considered: T , T a will denote the problem posed on a disk or a ring with (D) b.c. and T N , T aN , when posed with (N) b.c.Introducing the scaling r = t/R e and the function S(r) = T (rR e ), we write previous equation as According to the domain and the b.c., the solution will be denoted S, S N , S a and S aN as in the previous case.

Numerical issues
In this article we shall describe a Hybrid Symbolic-Numeric algorithm based on the Shooting Method, see Lemma 4.4, and discuss some qualitative properties exhibited by the numerical solutions of equation (3).One of these being the existence of thresholds, see Theorem 1 and the other being the estimates for the saturation lengths, see section 5.1.To this end, we notice that equation (3) has a singularity at the origin which means that we shall shoot from r = 1 and the evaluation runs backwards.We also notice that setting the target is also part of the algorithm.
On the other hand, in view of the saturation effect exhibited by the solutions, see section 5.1, solving the problem given by equation (3) appears to be very illconditioned.In Figure 1  We notice that, in any case, the absolute value for the error near the origin is equal to 1; another remarkable fact is that the divergence of solutions takes place in a very short interval, which is an alternative way to express the saturation effect.3) is ill-conditioned.

Setting the precision
Since parameters R i and R e enter into the equation ( 2) only on the boundary conditions and, as it is shown by Figure 1, for r 0 far from the boundary, the solution satisfies S(r 0 ) ∼ 1 and S (r 0 ) ∼ 0, it is natural to expect a strong relationship among solutions with different boundary conditions and with different values for the parameters R i and R e .In the sequel we shall exploit this fact to set the shooting parameter.Since we have based our algorithm on a high-precision arithmetic, we first explore the size of the working precision required in the computations and its relationship with the outer radius R = R e .In Figure 2 we plot the relative error for the pair of enclosing solutions produced in the first step of the algorithm as a function of R: Prec(R) = Log(ErrRel(R))/R, where ErrRel(R) = |S p (δ) − S n (δ)|S p (δ) −1 , see subsection 4.2 for details; in dotted line, the linear growth 0.7R.Motivated by these results, we set Prec(R) = 0.7R as the working precision.

The algorithm
Our algorithm is divided into two main steps.First, we shoot from r = 1 to some auxiliary target and advance by means of the Bisection Method to r = 0.Then, we adjust the solution to get the correct value using the Regula-Falsi Method.Since we are interested in four (highly related) problems (disk, ring, (D) and (N) b.c.) it is natural to expect (minor) differences when presenting them: each of one will be treated separately.We first discuss the problem posed on a disk and then we move to the problem posed on a ring.To handle the singularity at the origin, when dealing with the problem posed on a disk, we need to get some estimate for the solution for r ∼ 0.
Lemma 1 (Behaviour close to the origin) Let R > R c and let S(r) be a solution for (3) with S(0) = 0.Then, there exists a constant C(R) such that S(r) ∼ C(R)J 2 (Rr), valid for r ∼ 0.
Proof Let S(r) be the solution of (3).Since S(0) = 0 we have that S(r) ∼ 0, for r ∼ 0; thus, we have 1 − S(r) 2 ∼ 1.More precisely, let u be some rounding unit and let fl(•) be the related floating point representation, for |S(r)| u 1/2 we have fl(1 − S(r) 2 ) = 1 and the non linear term evaluates as the linear one: fl(S(r) 1 − S(r) 2 = fl(S(r)).This means that fl(S(r)) solves the linear equation: From previous lemma, for δ ∼ 0, the solution satisfies S(δ) ∼ C(R)J 2 (Rδ), in which C(R) is to be fitted: the second step of our algorithm deals with this issue.The first step consists in to produce two enveloping solutions that are defined in [δ, 1].This is achieved by shooting to the point (δ; 0): we look for solutions of (3), S p and S n such that S p (δ) > 0 and S n (δ) < 0. The main fact is that δ depends upon some estimate for the constant C(R).To track the behaviour in both limits and as R → R + c and as R → +∞, see Figure 3a, we express the constant as . This expression, together with the inequality J 2 (s) < s 2 /8, valid for s ∼ 0, allow us to deduce the estimate S(r) < (R − R c ) 1/2 R 3/2 r 2 /2, which we use to get uniform estimates in r and, therefore, to ensure the asymptotics S(r Clearly, δ = δ(R) and its actual value, according to Theorem 1, depends upon the problem to be considered: (D) or (N) b.c.We finally note that, since M (R) is to be found in the second step, we also need to set an auxiliary target for the fist step: we set S aux (δ) = 0.
Lemma 2 (Setting the target) Let 0 < u 1 be some arbitrary precision.For R > R c , and accordingly with the boundary conditions considered, we set δ and the target S(δ) as follows.Fig. 3 4.2 First step: Estimating the shooting parameter Since our algorithm proceeds by means of the bisection method, we need to have an upper and a lower solution, S p and S n respectively, that enclose the solution.To this end, we start considering equation ( 3) with Dirichlet boundary conditions.Let R > t 1 , we have the following empirical formula for the shooting parameter, in which In Figure 3b we show in red dots the actual value for S (1) and in solid blue line, the graph of |V (R)|.This estimate together with the monotonicity of )M α that yield the enclosing solutions needed to start the bisection method.
The algorithm proceeds by calling the Solver for equation (3) with S(1) = 0 and S (1) = v, with v = (v P + v N )/2 (bisection method) as the initial data; we notice that when calling the solver one must set Prec = 0.7R as the working precision: Solver has to have the option to set an arbitrary working precision.To enhance the performance, we force the solver to exit the computation as soon as the solution S reaches either S(dd) = 1 or S(dd) = 0; this gives an approximate solution defined in [dd; 1] which is updated as the upper solution S p or lower solution S n depending whether S(dd) = 1 or S(dd) = 0.The loop is designed to advance (backwards) from dd = 1 to dd = δ, with δ(u, R) given by Lemma 2, which acts as the cut-off criterion.The output is a vector containing what is needed to start the second step: {−v P , S p (δ), S (δ), −v N , S n (δ), S n (δ)}, the initial velocity, the computed solution and its derivative in δ for both the upper and lower solutions, respectively.We notice that, we also ask the solver for the value of S (dd) at each stage of the loop.

Second step: Adjusting the constant
We now move to the second step, which is to adjust the constant M (R).This is achieved by means of the Regula-Falsi method.With this in mind, we recall that the first step provides a solution S(r) for equation (3) which is defined in [δ; 1].We are looking for both an initial datum v R (the derivative at r = 1) and a constant M (R), such that the solution S(r) satisfies both M (R)J 2 (Rδ) = S(δ) and M (R)RJ 2 (Rδ) = S (δ) (in order to have a regular solution).This yields the condition S(δ)R = J 2 (Rδ)/J 2 (Rδ)S (δ), which is expressed as a function of the shooting parameter v: Thus, we use the Regula-Falsi Method to find a solution for the equation ( 4).We notice that the first step is designed in such a way that it produces outputs satisfying S P (δ) • S N (δ) < 0. Since S(dd + ) > 0 at each stage, the lower solution always satisfies S n (δ) > 0; however, since there is a small set of initial data for which the related solutions satisfy S(δ) > 0 and S (δ) > 0, we need to avoid this situation: this is accomplished introducing an auxiliary variable v 1 , that tracks the false upper solution; the algorithm proceeds by taking v = (v 1 + v N )/2, and updates the value in v P only when S(δ) > 0 and S (δ) < 0.

The pseudo-code
Below, we write the pseudo-code for each of the steps of the algorithm discussed above, that solves equation ( 3 . % The computed solution.

Behaviour of the solutions on the disk
The main result in this section is the presence of a saturation effect, where the radial distribution function converges to an almost constant solution close to 1 in a large interval.Conversely, the spatial variation of S(r) or T (r) in a (small) interval depends on the typical radii of the hollow cylinder and constitutes a purely geometrical effect.
If we poses the problem in a disk, this fluctuation also occurs around the origin since we require continuity when r → 0; when (D) b.c. are considered, same applies to the boundary point r = 1 or t = R (depending whether we are considering solutions of equation ( 3) or ( 2)), and we have the domain split into three regions.Finally, it should be noted that this effect becomes more visible as the radius becomes larger and is already present for R 100.See Figure 4.  We notice that, applying the scale t = Rr, the same results are valid to solutions for equation (2), posed in [0; R], the maximum being located at (Rx m (R); S m ).

Saturation effect
This subsection is devoted to discuss the saturation effect exhibited for the solutions of equation ( 3) posed on a disk with both (D) or (N) b.c.As indicated above, this feature allows us to divide the domain into two or three regions, each with qualitatively different behaviour, which will be considered separately.Finally, let us mention that we make use of this feature to get an approximate solution valid for R 1, compensating the fact that the algorithm becomes slower as R becomes greater.We first deal with the outer region.For this purpose we will consider equation ( 3) posed in [0; 1] (disk) with (D) b.c.Let R 1, let S be the solution of (3) and let also x m = x m (R) be the point in which S attains its maximum.The numerical evidence shown in Figure 4a suggests that the outer region [x m ; 1] shrinks as R 1. Thus, taking r ∼ 1 and setting g(r) = 1 − S(r + x m ), we get the auxiliary initial value problem, g (r) + g (r) + 4(1 − g(r)) = 2R 2 g(r)(1 − g(r))(1 − g(r)/2), with the initial data g (0) = −S (x m ) = 0 and g(0) = 1 − S(x m ) = g m that needs to be estimated.Since g (0) = (1 − g(0))(2R 2 g(0)(g(0)/2 − 4), we take g m ∼ g m := 2R −2 (1 + 2R −2 ), that leads to g (0) ∼ (1 − g(0))R 2 g(0) 2 > 0 which expresses the fact that r = 0 is a minimum.(Notice that, the simpler values g m = 0, 2R −2 are not useful.) Finally, we make the scaling t = Rr and obtain the approximate equation for p R (t) = (t/R): together with the initial conditions p R (0) = g m , p R (0) = 0, which can be solved numerically for values of R of large size, we tested it for R 5000 Thus, solving the equation p R (t) = 1 and calling s m (R) such that p R (s m (R)) = 1, we get both an estimate for x m in the form x m ∼ x m (R) = 1 − s m (R)/R and also an approximation for the solution in [x m ; 1] in the form: We show the relative error for the solutions Err   Regarding the problem to measure the saturation effect, that is to get an estimate for the length of the interval in which the solution S passes from S(1) = 0 to S(1 − η e ) ∼ 1, we denote by η e this saturation length and take as a (practical) definition the property S(1 − η e (R)) = 0. 9995 or, in terms of the auxiliary problem, p(Rη e (R)) = 0. 0005.Since we are also interested in provide an estimate for x m (R) for larger values of R, and since this is given by p(s m ) = 1, we treat both problems at the same time.
Remark 2 (Estimating the maximum and the (outer) saturation length ) Let (x m (R), S m (R)) be the location of the maximum for the solution S and let η e (R) denote the (outer) saturation length.Figure 6a shows the behaviour for R(1 − x m (R)) as a function of R, we notice the logarithmic growth; we also report the values obtained for the saturation length in Figure 6b.Below, we record the asymptotics: η e (R) = η * e R −1 with η * e ∼ 5. 8652.As we reported in Figure 4b, all the solutions of equation ( 2) with R 1 (for both (D) or (N) b.c.) tend to coincide when evaluated close to the origin; in addition, for R > 100 this kind of saturation effect is already present.As we have done before, we shall provide an estimate for the (inner) saturation effect and we also take profit of this feature in order to build an approximate solution for greater values of R.
Thus, we take R = 150 and denote by T 150 (t) the solution of equation ( 2 In Figure 7a, we plot the relative error (S(r)−T 150 (rR))/S(r) for r ∈ [0; 2η * 0 /240], where S is the computed solution for R = 240.The value η * i = 63.2693 is related with the saturation length of the solutions near the origin, discussed below.In addition, we note that the peak of order 10 −8 near r = δ reveals the value u = 10 −8 of the rounding unit used when setting the target, see Lemma 2.
When the problem is posed in a disk with (D) b.c., and look for bounded solutions, r = 0 acts as a boundary and it is worth to estimate the related saturation length; indeed, we are interested in make a comparison between these saturation lengths.
To estimate the (inner) saturation length presented in the solutions we look for η 0 (R) such that S(η 0 (R)) = 0. 9995.Thus, we take R = 150, and consider η * 0 = 63.2693 the smaller solution of the equation T 150 (t) = 0. 9995; notice that this gives an uniform length when the problem is posed in [0; R].Introducing the scaling t = rR, we obtain the saturation length for solutions of equation ( 3) as: η 0 (R) = η * 0 R −1 ; we also note that this value is more than ten times larger than η * e .
Finally, we notice that the values reported for the saturation lengths strongly depends on the choice 0.9995 involved in its definition, which in turn is completely arbitrary.However, in Figure 8a and 8b we plot the values related to the choice 0.995, which reveals some bias to the origin since η * e2 /η * e ∼ 1 while η * 02 /η * 0 ∼ 2.  A straightforward computation gives: h(r) = (x m R) −2 ((x m r) 2 + (x m r) −2 ), r = exp s, which could be simplified using the estimate x m (R) ∼ 1.We write the approximate solution, with domain [2η * 0 /R; x m (R)], (D) b.c., or [2η * 0 /R; 1], (N) b.c., as: In Figure 7b we show the relative error Err rel (r) = |S(r) − S app (r)|/S(r) for r ∈ [2η * i R −1 ; x m (R)], for R = 200, 240, 280; we notice that the error decreases as R increases.
Lemma 3 (Asymptotic solution on a disk with (D) b.c.)For R > 150 we get the following approximate solution, which is expressed using the functions introduced previously, we refer to Figures 5a, 7b, 7a for details of the relative errors involved.Here, we collect the minor changes that have to be made when dealing with Neumann boundary conditions on the disk.We recall that S, T shall denote a solution for equation ( 3) or ( 2) with (D) b.c. and S N , T N , a solution on a disk posed with (N) b.c.Since, both the working precision and the target have been set, the first of these changes has to do with the estimate of the shooting parameter and the setting of the enclosing solutions, for which we take profit of the strong relationship among the solutions shown in subsection 3.1 when posing the problem in [0; R]; however, it should be remarked that for R ∈ (s 1 ; t 1 ] there is no available solution for (D) problem, and only for R > 12, T serves as a guess, see the dotted red line in Figure 4d.Since we are considering Neumann boundary conditions, we take S N (1) ∼ S(x m (R)).To give a precise formulation, we recall the results of Theorem 1 concerning the threshold value, which now is given by R c = s 1 (the first positive zero of J 2 ), and also the estimate Remark 4 (Estimating the shooting parameter: (N) b.c.)Let R > s 1 , we have the estimate for the shooting parameter S N (1): Concerning the enclosing solutions, the result is similar to the one given previously.
For smaller values, say R ∈ [3.06; 12], try first with S p (1) = λ(R)W (R), interpolating between λ(3.06) = 1.6 and λ(12) = 1.05; on the other hand, the lower solution estimate S n = 0.99W (R). is still valid The first step proceeds with no further modification, except that when calling the solver one must set (N) b.c.The same applies to the second step.Therefore, the algorithm in this case is the one given in section 4.4 with the changes already mentioned.
Remark 6 (Asymptotic solution on a disk (N) b.c.)Let R 1 be fixed and let S N (r) be the solution of equation (3) with Neumann b.c.The approximate solution is given by: where S app (r) is given by Lemma 3 and x m (R) is the point in which S attains its maximum.Regarding the existence of thresholds for the inner radius R i , when the problem is posed in a ring the situation is as follows: for Dirichlet b.c. the problem has non trivial solutions for R i below a threshold R a (R e ), while for Neumann b.c. the problem has non trivial solutions for any small values of R i /R e , see Figures 9a and 9b.We discuss this issue in section 6, where we make a theoretical approach, see Theorem 1.The non-linear equation involved in the Regula-Falsi Method has to be reformulated, as it must now only reflect the boundary condition; we express it as a function of the shooting parameter v (we add a superscript to distinguish it from the solution in a disk):

Behaviour of the solutions on a ring
First, we consider the problem posed with (D) b.c.Following the ideas of section 5.1.1,mainly the solution p R of equation ( 5), we get an estimate for the outer region [x m (R e ); 1], when R e > 150.On the other hand, when R i > 150, and also R i R a (R e ) (far from threshold) same approach as before allows us to build an estimate for the inner region [R i /R e ; R i /R e (2−x m (R i ))].For R i ∼ R a (R e ) use the asymptotic expression given by Theorem 2. For R i < 150 and R e > 150, apply the estimate to the outer region and solve the problem with the algorithm in [R i ; 150].
Remark 7 (Asymptotic solution on a ring: (D) b.c.)Let R a (R e ) R i > 150, and let p R be the solution of equation (5).We have the following approximate solution: Since the asymptotic provided by Remark 7 is obtained from p R we deduce an estimate for both saturation lengths, each of which depends on the related radius.
Regarding the problem posed with Neumann b.c., as we have already mentioned, there is a non-trivial solution for any R i ∈ (0, R e ) see Lemma 1. Due to the saturation effect, for R 1 the solution S aN has two different regimes: , in which S N a (r) is almost constant (tiny curl above the solution for the disk in Figure 10), and r ∈ [R i /R e (1 + η N i ); 1], in which S N a (r) ∼ S N (r).We also note η N i η * e .
In this section we collect some theoretical results concerning the existence of thresholds.We recall that G(x, y) = S(Rr)e i2ϕ solves the Guinzburg-Landau (G-L) equation (10) on the ring D Ri,Re := {(x; y) : R i x 2 + y 2 R e }, if and only if, S(r) solves the equation ( 11) in [R i /R e ; 1]: Below, we present the existence of three kind of thresholds: two of them, establish a lower bound for the (outer) radius R e of the disk, the third one provides an upper bound for the thinness of the ring when the problem is posed with (D) b.c.In addition, we perform an asymptotic analysis around these critical values.
Theorem 1 (Existence of thresholds) Let t 1 and s 1 be the first positive zeroes for J 2 (x) and J 2 (x) respectively; let z 1 (R) be the largest solution in (0; 1) of J 2 (Rz)Y z (R)− J 2 (R)Y 2 (Rz).Then, the following statements hold true: -(D), on a disk: For R e t 1 , S ≡ 0 is the unique solution for equation (11)  Proof Let G D , G N be a solution for the (G-L) equation ( 10) in D R with Dirichlet or Neumann boundary condition respectively.We denote both with G.We also denote with R c any of the values t 1 or s 1 .
In order to establish the existence of positive solutions, we recall that the eigenvalues for ∇ 2 posed on D R = { x R}, with either (D) or (N) b.c., are −λ 2 k,m R −2 (here J k (λ k,m ) = 0 or J k (λ k,m ) = 0).Since in our problem we have G = S(r)e i2ϕ , we are looking for solutions belonging to the invariant subspace related to the eigenvalues −λ 2 2,m /R 2 (positive zeroes of J 2 or J 2 ).Multiplying (G-L) by G, integrating in D R , and using the integration by parts formula, we get: Rearranging terms and using the optimality of the first eigenvalue, we deduce: 0 and, hence G ≡ 0. (For (N) b.c., this is valid in the orthogonal complement of constant solutions.) When the domain is the ring D Ri,Re = {R i x R e }, we only have to notice that the eigenvalues are given by the positive zeroes of J 2 (xR i )Y 2 (xR e ) − J 2 (xR e )Y 2 (xR i ), and z 1 = R a (R e ) is the smallest value of R i for which this function does not vanish in (0; 1).Finally, we notice that z 1 (R e ) and 1 are consecutive zeros for the Bessel function, and it is well known the estimate 1 − z 1 (R e ) = π/R e , see [6].

Conclusions and further applications
In conclusion, the numerical analysis of the Ginzburg Landau equation shows the influence of the boundaries in the radial distribution of the nematic parameter S(r).
In fact, the radial distribution saturates to its equilibrium value S 0 = 1(T 0 = 1) in a region that depends on the size and shape of the cylindrical domain considered.From a physical point of view, this saturation effect arises in the capacity of the radial distribution S(r) (T (r)) to reach the saddle point in the nematic state driven by a purely geometrical effect, but normalized in units of the (thermal) coherence length of the system.Moreover, in the case of an annular domain we observe that the saturation is not necessarily symmetric with respect to the origin and depends on the ratio between the internal and external radii of the sample.This asymmetric behaviour remains if we consider a disk and we require a continuous solution of S(r) (T (r)) at the origin, due to the ambiguity of the molecular alignment at the center of the sample.

4. 1
Solving the non linear problem in a disk.The first remarkable fact when solving the problem on a disk with (D) or (N) b.c. is the existence of a threshold R c such that for R e R c the problem only has the trivial solution.Clearly, the actual value for R c depends on the b.c.considered, see Theorem 1 below for details.
(a) Behaviour of M (R), for Rc R 80.(b) Estimate for |V (R)|, the shooting parameter.

5. 1 . 1
Saturation near the outer boundary (for (D) b.c.) rel (r) = |S(r) − S app (r)|, r ∈ [x m (R); 1], for different values of R in Figure 5a, notice that as R becomes greater the error decreases and x m (R) → 1. Figure 5b shows the relative error for the estimate of the maximum Err rel (R) = |x m (R)− x m (R)|, for values R ∈ [50; 280], we also notice that the error decreases as R increases.

( a )
Behaviour close to the origin.(b) Behaviour close to r = 1.

5. 3
Solving the non linear problem in a ring.Now, we consider equation (3) posed on [R i /R e ; 1] for 0 < R i with both (D) and (N) b.c.

Remark 8 (
Comparing saturation lengths) Let R a (R e ) R i > 150 and let η i , η e be the inner and outer saturation lengths for solutions of equation (3) posed in [R i /R e ; 1] with (D) b.c.We have: with Dirichlet condition S(0) = S(1) = 0. -(N), on a disk: For R e s 1 , S N ≡ 0 is the unique solution for equation (11) with Neumann condition(S N ) (0) = (S N ) (1) = 0. -(D), on a ring: For R i z 1 = R a (R e), S a ≡ 0 is the unique solution for equation(11) with Dirichlet b.c.S(R i /R e ) = S(1) = 0.In addition, for all R e > t 1 , we have 1 − z 1 (R e ) ∼ π/R e , as a lower bound for the thinness of the ring.

Lemma 1 (
Existence of non-constant solutions in a ring: (N) b.c.)Let R i ∼ R e , then there exists S N , a non-constant solution for equation (3) posed in [R i R −1 e ; 1] with (N) b.c.
Store it as a lower solution.Else:RF P = RF a, v P = RF v. % Store it as an upper solution.End If tol = |RF a|.% update the tolerance.