Self-similar Singularity of a 1D Model for the 3D Axisymmetric Euler Equations

We investigate the self-similar singularity of a 1D model for the 3D axisymmetric Euler equations, which is motivated by a particular singularity formation scenario observed in numerical computation. We prove the existence of a discrete family of self-similar profiles for this model and analyze their far-field properties. The self-similar profiles we find agree with direct simulation of the model and seem to have some stability.


Introduction and Main Results
Whether the 3D Euler equations develop finite-time singularity is regarded as one of the most important open problems in mathematical fluid mechanics. Interested readers may consult the survey [7] and references therein for more historical background about this outstanding problem. In this paper we investigate the self-similar singularity of a 1D model for the 3D axisymmetric Euler equations, which is motivated by a particular singularity formation scenario observed in numerical computation [10]. It is hoped that this work may help to understand the singularity formation of the 3D Euler equations.
According to the numerical results reported in [10], the solutions to (1.1) develop selfsimilar singularity in the meridian plane for certain initial conditions with no flow boundary condition at r = 1. The solid boundary and special symmetry of u θ , ω θ and ψ θ along the axial direction seem to make the flow in the meridian plane remain hyperbolic near the singularity point and be responsible for the observed finite-time singularity. A 1D model which approximates the dynamics of the 3D axisymmetric Euler equations along the solid boundary of the periodic cylinder r = 1 has been proposed and investigated by Hou and Luo in [8]. The finite-time singularity of this model is proved very recently by Choi, Hou, Kiselev, Luo, Sverak and Yao in [4]. Motivated by this new singularity formation scenario, Kiselev and Sverak [9] construct an example of 2D Euler solutions in a setting similar to [10]. In this example, the gradient of vorticity is proved to exhibit double exponential growth in time, which is known to be the fastest possible rate of growth for the 2D Euler equations. This example provides further evidence that the new singularity formation scenario reported in [10] is an interesting candidate to investigate the 3D Euler singularity.
Inspired by the work of [8] and [9], Choi, Kiselev, and Yao proposed the following 1D model (we call it the CKY model for short) [5] on [0, 1]: x w(y, t) y dy, (1.2c) w(0, t) = 0, ρ(0, t) = 0, ∂ x ρ(0, t) = 0. (1.2d) This 1D model, whose finite-time singularity has been proved in [5], can be viewed as a simplified approximation to the 1D model proposed by Hou and Luo in [8]. Like the 1D model of Hou and Luo, the CKY model approximates the 3D axisymmetric Euler equations (1.1) on the boundary of the cylinder r = 1 with The positiveness of ρ x (x, t) near the origin creates a compressive flow, which is the driving force of the finite-time singularity, and we will use this fact in our construction in section 2. Our numerical simulation suggests that this 1D model develops finite-time singularity in a way similar to that of the 3D axisymmetric Euler equations on the boundary of the cylinder as reported in [10], and the singular solutions to this model also develop self-similar structure. We investigate the self-similar singularity of this CKY model in this paper.
Since the finite-time singularity of the CKY model takes place near the origin, we consider the following self-similar ansatz with ρ, W and U being the self-similar profiles, Plugging these self-similar ansatz into equations (1.2) and matching the exponents of (T − t) for each equation, we get (1.5) c w = −1, c u = c l − 1, c ρ = c l − 2.
In this paper we study the existence and properties of solutions to the self-similar equations. A key fact for the CKY model is that the Biot-Savart law (1.6c) in the self-similar equations can be rewritten as a local relation with a decay condition, We first ignore the decay condition (1.8b), then the self-similar equations with (1.6c) replaced by (1.8a) become a nonlinear ODE system with singular RHS at ξ = 0, i.e., the RHS does not satisfy the Lipschitz condition. We construct solutions to this ODE system near ξ = 0 using a power series method, which can naturally overcome the singularity of the RHS at ξ = 0. The power series are unique up to a rescaling parameter for a fixed leading order of ρ(ξ) at ξ = 0, and can be extended to the whole R + by solving the ODE system. Then we prove that the decay condition (1.8b) determines the scaling exponents, and there exist a discrete family of c l , corresponding to different leading orders of ρ(ξ), such that the decay condition (1.8b) holds for the self-similar profiles we construct . We prove this part with the assistance of numerical computation and rigorous error estimation. Given the decay condition (1.8b), we further analyze the far-field behavior of these self-similar profiles and prove that the profiles are analytic with respect to a transformed variable at ξ = +∞.
An interesting fact for this 1D model is that self-similar profiles (1.4) exist only for a discrete set of the scaling exponent c l , corresponding to different leading orders of ρ(x, t) at x = 0. The self-similar profiles we find agree with direct simulation of the model and seem to have some stability in the sense that for fixed leading order of ρ(x, 0), the singular solutions using different initial conditions converge to the same set of self-similar profiles.
The self-similar profiles we construct are non-conventional in the sense that the velocity does not decay to 0 at infinity but grows with certain fractional power, correspondingly, the velocity field at the singularity time is Hölder continuous. Such behavior is also observed in the numerical simulation of the 3D Euler equations in [8], which is very different from the Leray type of self-similar solutions of the 3D Euler equations, whose existence has been ruled out under certain decay assumptions on the self-similar profiles [2,1,3].
Our method of analysis is of interest by itself. The existence result replies on the use of a power series method to deal with the singularity of the self-similar equations at the origin, and some very subtle and relatively sharp estimates of the self-similar profiles. The same approach can be taken to analyze the self-similar singularity of Burgers equation and get results similar to those obtained in this paper. We are currently investigating the possibility of extending this method to study the singularity of the 2D Boussinesq system. Another novelty in our analysis is the use of numerical computation with rigorous error control, which is an important step in establishing the existence of self-similar solutions. Our strategy to rigorously control the numerical error, including the numerical error of the integration scheme for an ODE system and the roundoff error introduced due to floating point operation, is quite general and can be used for other purposes.
The rest of this paper is organized as follows. In section 2, we construct the local selfsimilar profiles using a power series method and extend them to the whole R + . In section 3, we prove that the decay condition in the Biot-Savart law determines the scaling exponents in the self-similar solutions. In section 4, we prove the existence of self-similar profiles for different leading orders of ρ(x, t) at the origin. In section 5, we analyze the far-field behavior of the self-similar profiles. In section 6, we present our numerical results.

Construction of the Near-field Solutions
In this section, we ignore the decay condition (1.8b), i.e. we use the local relation (1.8a) to replace the Biot-Savart law (1.6c) in the self-similar equations. We use a power series method to construct the self-similar profiles near ξ = 0, which we call the near-field solutions. Then we prove that the local solutions constructed in this way can be extended to whole R + .
We have the following Theorem Theorem 2.1. For fixed c l > 2, there exist a family of non-trivial analytic solutions to the self-similar equations (1.6) (with (1.6c) replaced by (1.8a)) near ξ = 0, corresponding to different leading orders of ρ(ξ) at ξ = 0, s, which is defined in (1.9).
Proof. Assume Based on the local relation in the Biot-Savart law (1.8a), we have Plugging (2.1) into (1.6) and matching the k-th (k ≥ 1) order term ξ k , we get If initially the leading order of ρ(x, 0) at x = 0 is s, which is greater than 2 (1.2d), then according to (1.2b), s will remain as the leading order of ρ(x, t) as long as the velocity field remains smooth. And as we discussed in Section 1, ρ (s) x (x, t) should be positive near x = 0 to produce finite-time singularity. So in the corresponding self-similar profile (2.1a), there is To make (2.2a) hold for 1 ≤ k ≤ s, we require Since ρ s = 0, we require To make (2.2b) hold for 1 ≤ k < s, we require Since c l > 2, and And to make (2.2b) hold for k = s, we require For k > s, to make (2.2) hold, the coefficients ρ k and U k should satisfy which means the power series (2.1) can be determined inductively.
We can achieve this by choosing u 0 r and ρ 0 /u 0 small enough to make the last two hold, and then choosing r large enough to make the first two hold. For example, let Then the choice of (2.12) would satisfy (2.10). And we will use induction to prove that for all k ≥ s, For k = s, (2.13) holds by (2.10). Assume now for s ≤ k < n, (2.13) holds. Then for k = n ≥ s + 1, based on (2.9a) we have (2.14) |ρ n | ≤ n−1 m=s |U m ||(n − m + 1)||ρ n−m+1 | (n − s)(c l /s − 2/s) .
Using the induction assumption and the fact that ∞ m=2 1 m 2 ≤ 1, we have where we have used the fact n ≥ s + 1 in the second inequality and (2.10) in the third inequality. Thus (2.13) holds for ρ n . Based on (2.9b), we have Using the induction assumption and the fact that ∞ m=2 1 m 2 ≤ 1, we get where we have used (2.10) and the fact that n ≥ 3, n 2 /(n − 1) 2 ≤ 9/4 in the last inequality.
So we just proved (2.13) by induction, which implies the power series (2.1) converge in some short interval [0, 1/r). This completes the proof of Theorem 2.1.
Remark 2.1. We require c l > 2 in Theorem 2.1. If c l = 2, there exist only trivial solutions to (1.6). If c l < 2, then c ρ < 0 according to (1.5), which means ρ(x, t) blows up in finite time according to (1.4). This is impossible since ρ(x, t) is transported by the fluid.

Remark 2.2.
For fixed c l , we have one degree of freedom ρ s > 0 in constructing the nearfield solutions (2.1), which can be easily verified to play the role of a scaling parameter (1.7). We will simply choose ρ s = 1 in our argument for the rest part of this paper.
The power series (2.1) we construct only converge in a short interval near ξ = 0. However, these local self-similar profiles can be extended to +∞.
Theorem 2.2. For c l > 2, the analytic solutions (2.1) we construct can be extended to the whole R + , resulting in solutions to the self-similar equations (1.6) with (1.6c) replaced by (1.8a). Moreover, we have for ξ > 0, Proof. Since c l + U 1 = (c l − 2)/s > 0, ρ s > 0, W s = (s − 1)U s > 0, based on the leading orders of the power series (2.1), we can choose ǫ < 1 r small enough such that Then we consider extending the self-similar profiles from ξ = ǫ to +∞ by solving an ODE system with initial conditions given by the power series (2.1). LetŨ (ξ) = c l ξ + U(ξ), then according to (1.6),Ũ (ξ), ρ(ξ) and U(ξ) satisfy the following ODE The right hand side of (2.20) is locally Lipschitz continuous forŨ(ξ) = 0, ξ = 0, so we can solve the ODE system from ǫ and get its solution on interval [ǫ, T ). We first prove that W (ξ) is positive on [ǫ, T ). Otherwise denote ξ = t as the first time W (ξ) reaches 0, i.e.

Determination of the Scaling Exponents
In our construction of self-similar profiles in the previous section, we did not consider the decay condition (1.8b). In this section, we prove that the decay condition determines the scaling exponent c l , i.e. only for certain c l does the decay condition hold.
Recall that for a fixed leading order of ρ(ξ), the self-similar profiles U(ξ), ρ(ξ) and W (ξ) depend on the scaling exponent c l only. So we can define a function G(c l ) as We will prove that G(c l ) < +∞ and it is a continuous function of c l . Then the existence of c l to satisfy the decay condition (1.8b) will follow from the Intermediate Value Theorem if we can show that there exist c l l and c r l such that Theorem 3.1. For a fixed leading order of ρ(ξ), s, and a scaling exponent c l > 2, construct the power series (2.1) with ρ s = 1 and extend the profiles to R + . Then the limit and G(c l ) is a continuous function of c l .
We first make the following change of variables, Then we have According to (2.5), (2.18) and the fact thatÛ (η) is monotone increasing, we have Before proving Theorem 3.1, we will first prove the following two Lemmas.
Proof. Assume that for some c l > 2, G(c l ) ≤ −2. Then according to (3.7) and the fact that U (η) is increasing, we have for all η > 0, Then we get It follows from (3.6a) that By direct integration and (3.7), we have for η large enough, Using this estimate in (3.6b), we get This implies Then we have for η large enough, Using this lower bound in (3.6c), we will get The constants C in the above estimates are positive and independent of η. The inequality (3.15) implies thatÛ (η) → +∞ as η → +∞, which contradicts with G(c l ) ≤ −2. This completes the proof of Lemma 3.1.
We add a subscript c l to indicate the dependence of the profiles on c l for the rest part of this section: Proof. We only need to prove that for fixed c 0 l > 2,Û c l (η),ρ c l (η) andŴ c l (η) as functions of c l are continuous at c l = c 0 l . In our construction of the power series using (2.9), we can easily see that the coefficients U k and ρ k depend continuously on c l . And based on the condition (2.10), there exist uniform upper bounds of these coefficients for c l in a neighbourhood of c 0 l . This means there exists a fixed ǫ small enough, such that W c l (ǫ),ρ c l (ǫ) andÛ c l (ǫ) are continuous at c 0 l . Then we use the continuous dependence of ODE solutions on initial conditions and parameter to complete the proof of this Lemma. Now we begin to prove Theorem 3.1. We use an iterative method which enables to get shaper estimates of the profiles after each iteration. We finally attain thatÛ c l (η) converges uniformly to G(c l ), with which we can complete the proof of this Theorem.
Proof. Consider c 0 l > 2, we will prove that G(c 0 l ) < +∞, and G(c l ) is continuous at c l = c 0 l . According to Lemma 3.1 and Lemma 3.2, there exists η 0 large enough and a neighborhood of c 0 Then for c l ∈ I 0 and η > η 0 , there exist ǫ 2 > 0, such that Using this in (3.6a), we have for c l ∈ I 0 and η > η 0 , Using direct integration and Lemma 3.2, we have for c l ∈ I 0 , η > η 0 , Using this upper bound ofρ(η) in (3.6b), we have for c l ∈ I 0 , η > η 0 , The first term in (3.22) is negative according to (3.7) and the second term is integrable for η > η 0 . Then using Lemma 3.2, we have for c l ∈ I 0 , η > η 0 , Putting this upper bound in (3.6c) and using Lemma 3.2, we get for c l ∈ I 0 , η > η 0 , Putting this upper bound ofÛ (η) back in (3.6b), we have for c l ∈ I 0 , η > η 0 which by direct integration gives that for c l ∈ I 0 , η > η 0 , Thus we have for c l ∈ I 0 and η > η 0 , Using this sharper upper bound of W (η) in (3.6c), we get for c l ∈ I 0 , η > η 0 , (3.28)Û c l (η) < C 9 ln ln η.
Using the upper bound ofŴ c l (η) (3.31) in (3.6c) we conclude thatÛ c l (η) converges uniformly as η → +∞ for c l ∈ I 0 and complete the proof of this Theorem.
To complete the proof of our main result Theorem 1.1, we still need to verify condition (3.2) for different s. And we leave this part to section 4.

Existence of Self-Similar Profiles
In this section, we verify that condition (3.2) holds for several s ≥ 2, i.e., there exist c l l and c r l > 2, such that G(c l l ) < 0, G(c r l ) > 0, with which we can complete the proof of Theorem 1.1. We will use the following Lemma, which allows us to prove (3.2) using only estimates of the profiles at some fixed η 0 . Lemma 4.1. Consider solving equations (3.6) with initial conditions given by power se- Proof. Since G(c l ) = lim η→+∞Û (η), andÛ(η) is increasing according to (3.6c) and (2.18), so if u 0 > 0, then G(c l ) > u 0 > 0, and we finish the first part of the Lemma (4.1b).
We will use numerical computation and rigorous error estimation to verify condition (4.1a) or (4.1c). We first numerically construct the power series (2.1) and then extend the local self-similar profiles to some η 0 by numerically solving an ODE system. Note that in equations (3.6), there are 1/η 3 terms in the right hand side, which can be very large for small η. This will make the numerical solutions sensitive to roundoff error. So instead of solving (3.6) directly, we make the following change of variables, (4.5)ρ(η) =ρ(η)η −5 ,W (η) =Ŵ (η)η −1 ,Ũ (η) =Û(η), and the equations satisfied by this new set of unknowns are given bỹ After this change of variables, there are only 1/η terms for this new system of ODEs.
We demonstrate how to rigorously bound the numerical error of the self-similar profiles at η 0 through the case s = 2, but the same procedure can be applied to other s to prove the existence of self-similar profiles. For the case s = 2, we choose c l l = 3 and c r l = 8 in (3.2).

4.1.
The case s = 2, c l = 3. We need to go though the following several steps.
Step 1 We need to bound the truncation error of the power series (2.1). To numerically compute the power series (2.1), we first truncate the power series to certain terms, and the truncation error can be bounded using (2.13). For the case s = 2, c l = 3, it can be easily verified that the following choice of ρ 0 , u 0 and r makes (2.10) hold: (4.7) u 0 = 1 9 × 162 , ρ 0 = 1 9 × 9 × 162 , r = 162.
Based on (4.5) and (2.1), at ξ = 10 −3 , corresponding to η s = 10 −1 , we have Using estimate (2.13), if we truncate the series (4.8) at k = 20, the truncation error for all the three series can be bounded by 10 −15 . We will numerically compute and use them as initial conditions to solve (4.6).
Step 2 We need to bound the roundoff error in computing the truncated power series (4.9). Denote f l as the floating point operation and assume a and b are two floating point numbers, which can be stored exactly on computer. Then by the IEEE standard rounding off rules [12], we have if f l(a ⊙ b) = 0, where ⊙ can be +, −, × and ÷, and ǫ is the machine precision. When using double precision floating point operation on matlab, there is (4.11) ǫ < 3 × 10 −15 .
The following two Lemmas will be intensively used to control the roundoff error.
The RHS of the inequaltities (4.12) only involve floating point operation, which implies that even though we cannot get the exact answer because of the roundoff error, we can get rigorous lower and upper bounds of the answer using only floating point operation. The basic idea in the above lower and upper bounds is compensating the roundoff error by multiplying (1 ± nǫ). Since ǫ is very small, these lower and upper bounds are actually very tight. We only prove the upper bound in (4.12a) for the case ⊙ is '+'. Others are similar.
Proof. According to the rounding rule (4.10), In the last inequality we have used |δ 1 |, |δ 2 |, |δ 3 |, |δ 4 | ≤ ǫ.  Let where ⊙ can be +, −, × and ÷. When ⊙ is ÷, Lemma 4.3 is obvious, so we omit its proof here. When computing the truncated power series (4.9), we need to first compute the coefficients U k and ρ k based on (2.9). The relation (2.9), Lemma 4.3 and Lemma 4.2 together allow us to get lower and upper bounds of the coefficients U k and ρ k inductively: Assuming now we have got lower and upper bounds of U m and ρ m for m < k, which are all floating point numbers, then based on (2.9a), and using Lemma 4.2 and Lemma 4.3 in each single step of the computation, we can get lower and upper bounds of ρ k . Then we use the same technique in (2.9b) to get lower and upper bounds of U k .
After getting the upper and lower bounds of ρ k and U k for all k ≤ 20, we use Lemma 4.3 and Lemma 4.2 in computing (4.9), and get lower and upper bounds of the truncated power series. Based on our computation, the difference of the lower bounds and upper bounds are all less than 5 × 10 −14 , which means if we use the upper bounds to approximate the real values of the truncated power series (4.9), the numerical errors are less than 5 × 10 −14 . We denote the upper bounds of the truncated power series (4.9) as (4.20)ỹ 0 = (w 0 ,ũ 0 ,ρ 0 ) T .
Since the truncation error of the power series are less than 10 −15 , we get a rigorous error bound of the self-similar profiles (4.20) at η s = 0.1, which we denote by
Based on the previous step,ỹ 0 and E 0 given by (4.20) and (4.21). Our approach to bound the error introduced in numerically solving (4.6) is simultaneously trackingỹ n and E n , i.e. in each step of the forward Euler method we updateỹ n and E n together, (4.25) (ỹ n , E n ) → (ỹ n+1 , E n+1 ).
So we get a rigorous error bound at the n + 1-st step, where I 1 = |I +Ah|E n is the propagation of error from the previous step, I 2 = h 2 /2|y ′′ (y i , x i )| is the local truncation error, and I 3 is the roundoff error in computingỹ n+1 in (4.26a).
Step 4 We first consider the first part I 1 . We can get lower and upper bounds of y 4 using |y 4 −ỹ n | ≤ E n , Lemma 4.3 and Lemma 4.2. Since A is explicitly given by (4.29), we can get lower and upper bounds for each entry of I + Ah using Lemma 4.3 and Lemma 4.2. The strategy is the same as Step 2 and we omit the details here. Taking absolute value for these lower and upper bounds, we get upper bounds for |I + Ah|. Then we use Lemma 4.3 and Lemma 4.2 in computing |I + Ah|E n and finally get rigorous upper bound of I 1 .
Step 5 Then we consider the second part I 2 . To bound the local truncation error I 2 , we need to bound the second order derivatives of the solutions on the interval [x n , x n+1 ]. For c l = 3 we haveρ We need the following a priori estimates to bound (4.33).
Lemma 4.4. Consider the ODE system (4.6) with initial conditions given by power series (2.1), assume at η 0 , the solutions areŨ(η 0 ),W (η 0 ),ρ(η 0 ), then for η ∈ [η 0 , η 0 + h], we have the following a priori estimates Proof. According to (4.6a) and the lower bound ofŨ(η) (3.7), we have By direct integration, we can get ρ max and ρ min . SinceŨ (η) is increasing according to (4.6c), so we get the lower bound u min . Then using the upper bound ρ max and the lower bound of U (η) (3.7) in (4.6b), we get By direct integration we get the upper bound w max . Putting the upper bound ofW (η) in (4.6c), we get the upper bound ofũ(η), u max . Using the upper bound w max and lower bound u min in (4.6b), we have and with this we can get the lower bound ofW (η), w min .
Using |y n −ỹ n | ≤ E n Lemma 4.3 and Lemma 4.2, we can get lower and upper bounds of y n . Putting these the upper and lower bounds in (4.34), and using the same strategy as Step 2, we can get lower and upper bounds of the solutions on the interval [x n , x n+1 ]. Finally using these lower and upper bounds of the solutions in (4.33) and using the same strategy as in Step 2, we can rigorously bound I 2 using only floating point operation.
Step 6 The last part we need to control is I 3 . I 3 is the roundoff error e roundoff in updating y n+1 using (4.26a). Using Lemma 4.3 and Lemma 4.2 and the same strategy as Step 2, we can get lower and upper bounds of the exact numerical solutions at the n + 1-th step using only floating point operation. In our computation we use the upper bounds as our numerical solutions at n + 1-th step, then the roundoff error introduced in this step can be bounded by the difference between the lower and upper bounds. So using Lemma 4.3 and Lemma 4.2 we can get upper bound of I 3 using only floating point operation.
Step 7 Putting the three parts of error bounds together and using Lemma 4.3 and Lemma 4.2, we get a rigorous error bound in the next step E n+1 using only floating point operation.
Step 8 Keep updatingỹ n and E n following Steps 3-7, we finally get the numerical solution of (4.6) at η 0 = 3 and a rigorous error bound, which are With these numerical solutions and the error bounds, we can easily verify that condition (4.1c) holds, and we complete the proof that for s = 2, The case s = 2, c l =8. The verification of G(8) > 0 can be done in the same way. In the construction of the local solutions, we can easily verify the choice of (4.41) u 0 = 1 6 , ρ 0 = 1 18 , r = 6, makes the constraint (2.10) hold. Then we truncate the power series (2.1) to the first 20 terms and evaluate them at ξ = 0.6 8 , corresponding to η s = 0.6. Using the same technique as the case c l = 3, we can get the approximated profiles and error bounds at η s , (4.42)ỹ 0 , E 0 = (10 −13 , 10 −13 , 10 −13 ) T .
Then we begin to solve (4.6). Using the same technique as the previous case to bound the numerical error in solving the ODE system, we finally get With G(3) < 0, G(8) > 0, we conclude that there exists a c l such that the self-similar equations have solutions, and the leading order of ρ(ξ) at ξ = 0 is s = 2.
Following the same procedure we also verified that for Remark 4.2. We only verify the existence of self-similar profiles for s = 2, 3, 4, 5. But the same procedure can be applied to s > 5 to verify the existence of self-similar profiles.

Behavior of the Self-Similar Profiles at Infinity
In this section, we prove that if the decay condition (1.8b) in the Biot-Savart law is satisfied for the self-similar profiles we construct, then the profiles are actually analytic with respect to a transformed variable θ = ξ −1/c l at θ = 0. With this we can complete the proof of Theorem 1.2. This far-field property of the profiles can explain the Hölder continuity of the velocity field at the singularity time observed in numerical simulation of this model. Theorem 5.1. For some fixed c l and s, if the self-similar profiles constructed using power series (2.1) and extended to whole R + satisfy the decay condition (1.8b), then after the following change of variables, U (θ),W (θ) andρ(θ) are analytic functions at θ = 0.
Our strategy is the following: We first prove that after the change of variables (5.1),Ũ (θ), W (θ) andρ(θ) are smooth at θ = 0. Then we argue that there exist analytic solutions to the ODE system they satisfy with the same initial conditions at θ = 0. Finally we prove that smooth solutions with given initial conditions are unique and complete the proof.
Using a simple bootstrap argument, we can getW (θ),ρ(θ) andŨ (θ) are in C ∞ [0, +∞) . On the other hand, given the initial conditions (5.12d), we can construct the following power series solutions to equations (5.12): Plugging these power series ansatz in (5.12) and matching the coefficients of θ k , we can uniquely determine the coefficientsŨ k ,W k ,ρ k and prove that the power series (5.14) converge in a small neighborhood of θ = 0. We omit the details here, because the argument are the same as in section 2. Then to prove the analyticity ofŨ(θ),W (θ) andρ(θ) at θ = 0, we only need the uniqueness of smooth solutions to (5.12) with initial condition (5.12d). The RHS of (5.12c) is not Lipschitz continuous, so the classical result will not apply here.
The C in the above estimates are all positive constants independent of ǫ. Choosing ǫ small enough, we get a contradiction in (5.19) and (5.20), thus which means the solution is unique. And we complete the proof of this Theorem. The above Theorem implies the self-similar profiles we construct are non-conventional in the sense that the velocity does not decay to 0 at +∞ but grows with certain fractional power. Coming back to the self-similar ansatz (1.4), we have For t close to T , based on Theorem 5.1 , we have This can explain the Hölder continuity of the velocity field near singularity time, which is observed in simulation of the 1D model. Such behavior is also observed in the numerical computation of the 3D Euler equations in [8], which is very different from the Leray type of self-similar solutions of the 3D Euler equations, whose existence has been ruled out under certain decay assumptions on the self-similar profiles [2, 1, 3].

Numerical Results
In this section we numerically locate c l which makes G(c l ) = 0 for several different s and construct the corresponding self-similar profiles. The scaling exponents and self-similar profiles obtained from solving the self-similar equations (1.6) agree with those obtained from direct numerical simulation of the CKY model. We also find that for fixed s, the singular solutions using different initial conditions converge to the same self-similar profiles after rescaling, which means the self-similar profiles we find have some stability property.
6.1. Numerical methods for solving the self-similar equations. For any fixed c l > 2, we first numerically compute the coefficients ρ k , U k in (2.1) up to k = 50 and determine the convergence radius of the power series using the following linear regression for s ≤ k ≤ 50, (6.1) log ρ k = k log r 1 + c 1 , log U k = k log r 2 + c 2 .
We choose and construct the local self-similar profiles on [0, r/2].
Then we continue solving equation (1.6) from ξ = r/2 to ξ = 1 using the 4th order explicit Runge-Kutta method with step-size h = 1−r/2 10 4 . After ξ = 1, we make the change of variables (3.4) and begin to solve (3.6) from η = 1 to η = 10 5 using 4th order Runge-Kutta method with step-size h = 10 5 −1 10 6 . We useÛ c l (10 5 ) as an approximation to G(c l ). We use the bisection method to find the root of G(c l ), and the stopping criterion is that the length of the subinterval becomes less than 10 −5 . After we get the scaling exponent c l , we construct the local self-similar profiles using power series (2.1) as before and numerically extend them from ξ = r/2 to ξ = 10 using the explicit 4th order Runge-Kutta method with step-size h = 9 10 4 . Then we locate the maxima of W , which is W max = W (ξ 0 ). We consider s = 2, 3, 4, 5, and for these cases ξ 0 are all less than 10. Finally we rescale the maximum of the W (ξ) to (1, 1), and get the rescaled self-similar profile of w We will only compare the self-similar profiles of W s with direct simulation of the CKY model in this paper, but the numerical results for the profiles of ρ and U are similar. 6.2. Numerical methods for simulating the model. In the direct simulation of the CKY model, we use a particle method. We consider N + 1 particles with position, density and vorticity given by In computing the velocity field, we use the trapezoidal rule to approximate (1.6c), In computing the driving force of w, which is ρ x , we use the three points rule: Initially, 10 5 +1 particles are equally put in the short interval [0, 10 −3 ], which are sufficient to resolve the solutions in the self-similar regime. Outside this short interval 10 5 − 10 2 particles are equally placed with distance 10 −5 . So the total number of particles is N = 2 × 10 5 − 10 2 .
Then we need to solve the following ODE system The initial condition of ρ is (6.8) ρ(x, 0) = (1 − cos(πx)) s/2 , whose leading order at x = 0 is s.
We solve the ODE system 6.7 using the 4-th order explicit Runge-Kutta method, and the time step dt is chosen adaptively to avoid the particles cross each other: Simulation of this model will stop once the maximal vorticity reaches some preset W max . And at each time step, we record the maximal vorticity w max (t i ), and the position where it is attained q max (t i ). According to the self-similar ansatz, we have Thus we can compute c l c w , and the singularity time T by doing linear regression, We compute the time derivatives of log w max (t) and log q max (t) using the center difference method, and the linear regressions are done in some time interval close to the singularity time while the numerical solutions still have good accuracy.
At certain time steps close to the singularity time, t i , i = 1, 2, 3, let w i be the maximal vorticity at time t i and q i be the position the maximal vorticity is attained. We rescale the numerical solution and get the self-similar profile of w, In the next subsection we will compare the self-similar profiles W i s (ξ) (6.12), which are obtained from direct simulation of the model, with W s (ξ) given by (6.3), which is obtained from solving the self-similar equations (1.6).
Near the singularity time the velocity field seems to be Hölder continuous near the origin, Then we can determine the Hölder exponent α by doing linear regression (6.14) ln u(x, T ) ≈ ln C + α ln x.
We will compare the exponents α we get from direct simulation of the CKY model with our prediction (5.23) obtained from analyzing the self-similar equations.
6.3. Comparison results. In our direct simulation of the CKY model, we first choose the initial condition of w(x, t) as We compute the scaling exponents c w and c l for different leading orders of ρ, s = 2, 3, 4, 5 using direct simulation of the CKY model, and the results are listed in Table 1 and Table 2.
We also compute the Hölder exponents of the velocity field near the origin at the singularity time and compare them with 1−1/c l as predicted by (5.23). The results are listed in Table 3. The c l we use are those obtained from solving the self-similar equations.
For the case s = 2, the dependence of G(c l ) on c l is plotted in Figure (1). From this Figure, we can see that G(c l ) depends continuously on c l as we have proved. Besides, G(c l ) seems to be a monotone increasing function, which implies that for fixed s, the scaling exponent c l to make the decay condition (1.8b) hold is unique. .  We also compare the self-similar profiles obtained from solving (1.6) and those obtained from direct simulation of the model. For different s, they are plotted in Figure 2. The lines labeled 'exact' are obtained from solving the self-similar equation. Others are obtained from rescaling the solution at different time steps corresponding to different maximal vorticity.
To demonstrate the stability the self-similar profiles, we consider another initial condition, (6.20) w(x, 0) = x − x 2 .
Again we compare the self-similar profiles of w from direct simulation of the model with those obtained from solving the self-similar equations (1.6). They are plotted in Figure 3.   From Figure 2, 3, we can see that after rescaling, the singular solutions at different time steps before the singularity time are very close, which implies that the solutions of this 1D model develop self-similar singularity. And the self-similar profiles obtained from direct simulation of the model (6.12) agree very well with the self-similar profiles we construct (6.3) by sovling the self-similar equations (1.6). For fixed leading order of ρ(x, 0) at the origin, the singular solutions with different initial conditions converge to the same set of self-similar profiles, which implies that the profiles we construct have some stability property.