1. Introduction
The basic mechanism of the thin-layer, shear-flow instability, known as the Kelvin–Helmholtz (KH) instability, refers to the fluid motion generated at a shear zone separating two uniform layers, of generally different density, moving in equal and opposite directions. If the shear layer is taken to have zero thickness, the flow is termed a ‘vortex-sheet’ instability. Helmholtz (Reference Helmholtz1868) first discussed the fluid dynamical prototype while Lord Kelvin (Thompson Reference Thompson1871) later developed the linear theory for vortex-sheet KH instability. The KH instability is important and often dominant in many atmospheric and oceanic fluid phenomena that often involve the transition to and production of turbulence.
Vortex-sheet instability is known to exhibit ill posedness in general, in the sense of unbounded growth of arbitrarily small-wavelength disturbances. Consequently, the fluid dynamics beyond the linear growth stage has been widely studied. See; Moore (Reference Moore1978), Krasny (Reference Krasny1986b), Baker & Shelley (Reference Baker and Shelley1990), Sohn (Reference Sohn2016) and Caflisch, Lombardo & Sammartino (Reference Caflisch, Lombardo and Sammartino2020). Rosenhead (Reference Rosenhead1931) performed a numerical simulation using twelve point vortices per wavelength, observing a tendency of the vortex sheet to roll up into vortex cores. Calculations with up to twenty vortices were discussed by Birkhoff (Reference Birkhoff1962), who reported the onset of non-smooth, random-like motion, behaviour later associated with loss of analyticity (Ely & Baker Reference Ely and Baker1994) in the vortex-sheet evolution.
An advance in understanding was made by Moore (Reference Moore1979), who showed that a small-amplitude initial disturbance to an initially flat vortex sheet whose evolution was governed by the Birkhoff–Rott (BR) equation (Rott Reference Rott1956; Birkhoff Reference Birkhoff1962) produced a curvature singularity in a finite (critical) time $t_c$ proportional to the logarithm of the inverse disturbance amplitude. Verification of this prediction has been the subject of several high fidelity numerical simulations (Krasny Reference Krasny1986a; Shelley Reference Shelley1992; Nitsche Reference Nitsche2001). Meiron, Baker & Orszag (Reference Meiron, Baker and Orszag1982) solved the BR equation using a high-order Taylor series, finding singularity formation in a finite time. Existence of solutions up to a near-singularity formation time for a slightly perturbed, flat vortex sheet was demonstrated by Caflisch & Orellana (Reference Caflisch and Orellana1986). Exact solutions of the BR equation that show singularity formation have been reported by Duchon & Robert (Reference Duchon and Robert1988) and also Caflisch & Orellana (Reference Caflisch and Orellana1989). Cowley, Baker & Tanveer (Reference Cowley, Baker and Tanveer1999) explored the dynamics of the singularity formation event, showing that it can be interpreted as the convergence of a $3/2$-power singularity in the extended complex plane of the Lagrangian vortex sheet parameter, onto the real axis at $t=t_c$. They also estimated the vortex-sheet shape as $t\to t_c^-$ from below, giving a local form of the curvature singularity formed by a perturbation from a locally flat sheet.
Vortex-sheet motion for $t>t_c$ has been studied using the vortex-blob formulation (Krasny Reference Krasny1986a; Baker & Pham Reference Baker and Pham2006; Sohn Reference Sohn2016), which removes singularity formation, allowing long-time numerical simulation. Blob-based simulations show fine-scale geometry near the singularity formation region for $t>t_c$ at scales below the blob size. Whether the sub-blob evolution scales faithfully depict the immediate post-singularity vortex-sheet evolution very near the point of singularity formation, or are a blob-size-dependent artefact of the blob approximation, remains an open question. See Caflisch et al. (Reference Caflisch, Gargano, Sammartino and Sciacca2022) for discussion and Baker & Pham (Reference Baker and Pham2006) for detailed comparative numerical simulations.
Alternative approaches study the limit when a small parameter approaches zero, of initial-value solutions of either the Euler or Navier–Stokes equations. Baker & Shelley (Reference Baker and Shelley1990) solve the Euler equations describing the motion of a thin, uniform vorticity layer whose thickness is small compared with the local layer-centreline curvature, observing the formation of a central bulge in the region where singularity formation is expected for the vortex-sheet limit. Within the bulge exists a double-branched spiral curve whose material initial condition corresponds to the initial centreline of the perturbed layer. Using numerical Navier–Stokes simulations, Caflisch et al. (Reference Caflisch, Gargano, Sammartino and Sciacca2022) report the formation of an inner core with internal properties scaling on the Reynolds number. Their intermediate Reynolds number simulations reveal the formation of dual, separated vorticity cores that later merge.
Both the blob and physical-parameter simulations indicate that, prior to the singularity formation event, there is convergence to the vortex-sheet limit when the regularization parameter approaches zero; see Baker & Pham (Reference Baker and Pham2006) for blob calculations. Finite-time analyticity of the BR vortex-sheet problem with analytic initial data has been established; see Sulem et al. (Reference Sulem, Sulem, Bardos and Frisch1981), Majda & Bertozzi (Reference Majda and Bertozzi1992) and Wu (Reference Wu2006). Owing to the singular nature of the initial condition generated by singularity formation, the problem of continuing the solution of the BR equation for $t>t_c$ may not be unique, resulting in the existence of several, many and possibly infinitely many admissible classical continuation solutions. Hence, convergence of regularized solutions to one or more classical BR solutions in the immediate vicinity of the singularity formation region remains uncertain, and may be regularization dependent.
Presently, we describe the construction of a solution of the BR equation as a viable continuation of the Moore solution, valid for arbitrarily small time beyond $t_c$. The arguments and analysis comprise Taylor-series expansion, series summation, analytic continuation of summed expressions and numerical solution of both nonlinear and linear forms of the BR equation in several different but overlapping solution regimes. A summary of the role played by these elements is as follows.
(i) The full solution domain consists of a single spatial period of an infinite, flat vortex sheet in an incompressible, inviscid fluid of uniform density, that is subject to a generalized, single mode, periodic disturbance of small amplitude $\epsilon$ at $t=0$. Defining $\tau =t-t_c$, then for $\tau >0$ this is partitioned into three sub domains:
(a) An outer domain denoted region I comprising the single-period sheet.
(b) An intermediate domain denoted region II. This is defined as the region of validity of a solution comprising the most singular term in the time-wise Taylor-series expansion of the BR equation, using the inner expansion of the Moore solution at $t=t_c$ as initial condition.
(c) An inner domain, denoted region III, that is close to, and includes the singular point generated at $t=t_c$.
Region III exists only for $t>t_c$. A proper description of region II requires a prior derivation of an appropriate solution. Quantitative definitions of the three sub-domains will be given at the end of § 2.
(ii) In § 2 we extend the series approach of Moore (Reference Moore1979) using a generalized, single-mode initial condition. The series can be summed to provide a small-$\epsilon$ solution for the full-sheet structure near $t=t_c$ in terms of polylogarithm functions. This enables determination of the leading-order (in the cumulative circulation variable $\varGamma$) sheet shape near the point of curvature singularity formation at $t=t_c$. This is termed an intermediate initial condition for a region-II solution, comprising a flat sheet with an added $\varGamma ^{3/2}$ term that contains the singular behaviour. The polylog solution can be analytically continued for $t>t_c$ outside the range of convergence of the generating Fourier series. This shows sheet-branch separation or tearing at the sheet symmetry point at $t=t_c^+$.
(iii) Next, in § 3, the first two terms of a Taylor series in time are found based on repeated time-wise differentiation of the BR equation for $t>t_c$. These give the sheet velocity and acceleration at $t=t_c$ that are compared in § 4 with numerical evaluation of the Biot–Savart induced velocity from the Moore solution. In § 5, this analysis is extended to give the full Taylor-series solution for the most singular term of the BR equation. This series is shown to have a finite radius of convergence in a similarity parameter $\eta \sim \varGamma /(t-t_c)$, but can be summed to provide a closed-form expression defined as a region-II solution. The radius of convergence is taken as the inner bound of region II. This solution can also be analytically continued inside the series radius of convergence towards $\eta \to 0$, which is equivalent to $\varGamma \to 0$ for $t>t_c$. Sheet tearing is also found, which is shown in an appendix to be in quantitative agreement with the analytic continuation of the extended Moore solution in region I.
(iv) An inner solution in region III is developed in § 6. It is noted that the analytically continued region-I and region-II solutions display free vortex-sheet endpoints or tips that, without endpoint vortex-sheet roll up, are considered non-physical. The tearing feature, together with the perturbed similarity form of the region-II solution, nonetheless suggest an inner similarity expansion of the BR equation, leading to a nonlinear equation that describes a dominant or zeroth-order inner solution, together with a first-order correction in the small parameter $(t-t_c)^{1/2}$, defined by a linearized BR equation. It is argued that a solution of the zeroth-order equation based on double-spiral, vortex-sheet roll up about the singular point does not exist. The only available alternative is a sheet geometry comprising a pair of spatially separated, antisymmetric algebraic spirals. This structure is interpreted as consistent with the sheet tearing suggested by the analytically continued region-II solution, while restoring the physically expected, end-sheet roll up.
(v) Numerical solutions to both the nonlinear zeroth- and linear first-order equations in region III are obtained that effectively match the corresponding region-II solution components for large $\eta$. When the similarity plane is transformed back to the physical plane using a compression factor proportional to $\tau$, the zeroth-order, spiral vortices show roll-up centres that mutually recede as $t-t_c$, together with a $(t-t_c)^{3/2}$ correction provided by the outer flow. Near the sheet endpoints, the order-$(t-t_c)^{3/2}$ solution component agrees qualitatively but not quantitatively with the analytically continued region-II solution.
(vi) The composite region-II–region-III solutions are discussed in § 6.4 while conclusions are summarized in § 7.
2. Moore solution in region I
2.1. Problem definition
First, we extend Moore's solution for the evolution of a spatially periodic vortex sheet with a family of initial perturbations. Evolution of the vortex sheet in two-dimensional potential flow is described by a parameterized complex curve
where $t$ is time and $\varGamma$ measures the circulation in the sheet between an arbitrary point with coordinate $z$ and a reference fluid particle chosen as $\varGamma =0$. The sheet shape evolution in time is governed by the BR equation
where ${\bar {z}}$ is the complex conjugate of $z$ and the right-hand side is understood as a Cauchy principal-value (CPV) integral. We will assume that $z(\varGamma,t)$ satisfies two conditions
Using the second condition, (2.2) can be written in the form
We consider an initial-value problem that uses a general initial sinusoidal perturbation of the form
with initial amplitude $0<\epsilon \ll 1$ and phase $\phi \in [0,2{\rm \pi} )$. The real component $\epsilon \cos \phi$ stretches the circulation distribution while the imaginary component $\epsilon \sin \phi$ perturbs the sheet shape. Moore (Reference Moore1979) studied the particular case where $\phi ={\rm \pi} /2$. See Krasny (Reference Krasny1986b), Cowley et al. (Reference Cowley, Baker and Tanveer1999) and Shelley (Reference Shelley1992) for discussion of and results for alternative initial conditions. The unperturbed sheet with $\epsilon =0$ corresponds to the uniform shear flow with velocity field given by
2.2. Fourier expansion
For $t>0$, Moore (Reference Moore1979) expresses the sheet function in a Fourier series
here with initial conditions according to (2.5)
It immediately follows from symmetry that
Substituting (2.7) into (2.2) leads to a system of ordinary differential equations (ODEs) for the Fourier coefficients $\mathcal {A}_n$
where $I_K$ is the principal-value integral with respect to $\xi =\hat {\varGamma }-\varGamma$
where $r_p=1,2,3,\ldots$ are indices of $\mathcal {A}$, and $p=1,2,3,\ldots$ are sub-indices for $r$.
To estimate $\mathcal {A}_n$, the system of (2.10) can be further partitioned by the power series expansion
giving subsystems of ODEs for each $\mathcal {A}_{n,0}$, $\mathcal {A}_{n,2}$, $\mathcal {A}_{n,4}$,…. In particular, at leading order, $O(\epsilon ^{|n|})$, the corresponding integral defined in (2.11) can be calculated using residues as
leading to the following evolution equations for $\mathcal {A}_{n,0}$:
The first two equations read explicitly as
2.2.1. Asymptotic coefficients
The ODE system (2.14) and (2.8) can be solved recursively to give, for example, the asymptotic ${\mathcal {A}}_{1,0}$ and ${\mathcal {A}}_{2,0}$ for large $t$ as follows:
where
Here, we require $\lambda _1\neq 0$, i.e. $\phi \neq {\rm \pi}/4, 5{\rm \pi} /4$, to retain the leading-order behaviour. Special cases with $\lambda _1=0$ are not considered in this study.
Moore (Reference Moore1979) showed by induction that the solution for the leading-order $\mathcal {A}_{n,0}$ takes the form
where, by abbreviating $\lambda _{n} \equiv \lambda _{n,0}$, the leading-order coefficient $\lambda _n$ for $n\geq 2$ is given by the following recursion:
To determine the large $n$ behaviour of $\lambda _n$, Moore (Reference Moore1979) proposed a generating function $g(x)$ defined for $x\in \mathbb {R}$ as
with
given by (2.17). The recursion (2.19) thus implies
Equations (2.22) and (2.21) can be solved in closed form as
where $W_0$ is the Lambert W function. Taylor expanding (2.23) around $x=0$ yields
Therefore, comparing (2.20) and (2.24) shows
Here, $c_n>0$ given in (2.24) also corresponds to the special case $\lambda _1=1$ (i.e. $\phi ={\rm \pi} /2$) studied by Moore (Reference Moore1979), where it is shown that
as $n\to \infty$. Finally, substituting (2.26), (2.25) and (2.18) into (2.12) completes the following asymptotic Fourier coefficients $\mathcal {A}_n$ for $t\gg 1$, $n\gg 1$ and $t\gg n$:
Compared with the special case of $\lambda _1=1$, it is seen that a general $\lambda _1\neq 0$ simply rescales the initial perturbation size $\epsilon$ to $|\lambda _1|\epsilon$ for the coefficients $|\mathcal {A}_n(t)|$.
2.2.2. Critical time of singularity formation
Owing to correction terms contained in (2.18), the leading-order $\mathcal {A}_{n,0}$ given by (2.27), as an approximation of $\mathcal {A}_{n}$, deteriorates in accuracy as $n\to \infty$ for fixed $t$. Moore (Reference Moore1979) resolved this non-uniformity by introducing a strained time $s$, and consequently obtained
where $\lambda _{1}=1$ and $\alpha _1=2$. In Appendix A, we show that (2.28) also holds for the general case of $\lambda _{1}\neq 1$, where the constant $\alpha _1=\alpha _1(\phi )$ as a function of $\phi$ is given by (A13).
Equations (2.12) and (2.28) imply that, in contrast to (2.27), the uniformly convergent Fourier coefficients are expressed in terms of $s$ as follows:
Therefore, at a critical time $t_c(\epsilon )$ defined by
the $\mathcal {A}_n$ values for the Fourier series (2.7) given by (2.29) lose their exponential decay for $t\ge t_c$, and transition to algebraic decay as $O(n^{-5/2})$. Consequently, the initially smooth vortex sheet spontaneously loses analyticity of its shape function $z$ at $t_c$, when a singularity develops. This singularity is of the same form as that exhibited in (2.27), but occurs at a slightly different time. For small perturbation size $\epsilon$, and consequently large critical time, $t_c$ can be well approximated by solving (2.30) at leading order to give
where $W_0$ is again the Lambert W function.
2.3. Sheet structure at critical time
Using (2.7), (2.8) and (2.12), the series describing the vortex-sheet shape can be asymptotically summed for small $\epsilon$ as
With $\mathcal {A}_{n,0}$ given by (2.28), the Fourier series (2.32) converges uniformly to
for all $t\leq t_c$ up to the critical time $t_c$ given by (2.31), where $\textrm {Li}_n(s)\equiv \sum _{p=1}^\infty s^p/p^n$ is the polylogarithm function. Evaluating (2.33) at $t=t_c$ yields the critical sheet shape, ${z_c=z(\varGamma,t_c)}$, as follows:
Equation (2.28) is valid for $n$ large, and the addition of low-mode corrections would modify (2.34). This is expected to be analytic and so will not change the nature of the singularity.
If $\lambda _1>0$ then, near the origin, where $0 \le \varGamma \ll 1$, (2.34) takes the asymptotic form
where
and $\zeta (3/2)\approx 2.61238$ is a constant following the Riemann zeta function $\zeta ({\cdot })$. Both (2.33) and (2.34) satisfy $z(-\varGamma,\tau ) = -z(\varGamma,\tau )$ as required but (2.35) does not because the complex constants are tailored to $\varGamma \ge 0^+$. This will not affect the subsequent analysis. An expression for the sheet shape at $t=t_c$, $\varGamma <0$ which satisfies the proper symmetry can be obtained but is not needed. Equation (2.35) agrees with a similar result at $t=t_c$ obtained by Cowley et al. (Reference Cowley, Baker and Tanveer1999), although estimates of constants differ.
Similarly, if $\lambda _1<0$, (2.34) can be expanded around $\varGamma ={\rm \pi}$ to show the same power series as (2.35). Since the role of $\lambda _1\neq 1$ in both cases is to rescale $\epsilon$ while preserving the same sheet shape near the singularity, we subsequently commit to the Moore (Reference Moore1979) initial value with $\lambda _1=1$, i.e. $\phi ={\rm \pi} /2$, without loss of generality. The variation of both the real and imaginary parts of $b(\epsilon )$ and $A(\epsilon )$ ($\lambda _1=1$) are displayed in figure 1. For $\epsilon \le 0.0439064$, the argument of $A$ lies in the fourth quadrant, crossing to the third quadrant for larger values. Later the special case $\mathrm {Re}[A]=\mathrm {Im}[A]$ or $\arg [A] = 5{\rm \pi} /4$ will arise. This occurs with $\epsilon = 0.248985,A = -0.639680 - 0.639680\textrm {i}, b=i$ which is nominally out of range of the present discussion. Calculations will generally use $\epsilon = 10^{-3}$ for which $b = 0.791266 + 0.208734\textrm {i}$, $A = 0.116149 - 0.199386\textrm {i}$. Asymptotic results with $\epsilon \to 0$ will also be presented.
Figure 2 shows the sheet shape in $0\le x \le {\rm \pi}$ for $\epsilon = 10^{-3}$ given by both (2.34) and the leading-order approximation (2.35). With $\epsilon = 10^{-3}$, the sheet angle at the origin at $t=t_c$ is $\arg [b(10^{-3})] = 0.257922$ radians or $14.7779^\circ$. The vortex-sheet strength, $\gamma \equiv |\partial z/\partial \varGamma |^{-1}$, which measures the jump in tangential velocity across the sheet, can be calculated from (2.34) at $t=t_c$ as
which is continuous for all $\varGamma \in [0,2{\rm \pi} ]$. From (2.37) or using (2.35)–(2.36a,b) it can be shown that, with $t=t_c$, as $\epsilon$ increases and $A$ crosses the imaginary axis at $\epsilon = 0.0439064$, then $\gamma _c$ changes from a local maximum, with negative singular gradient at $\varGamma \to 0^+$, to a local minimum while the sheet angle at the singular point passes through $\arg [b] = 45^\circ$. Since $\epsilon <<1$ for the Moore solution we will generally consider $\mathrm {Re}[A]>0$.
Using (2.35), the sheet curvature $\kappa$ for small $\varGamma$ follows as
with $\alpha =\sqrt {2/{\rm \pi} }\zeta (3/2)\approx 2.084$, which is singular when $\varGamma \to 0^+$.
2.4. Continued solution beyond critical time
We now address the motion of the sheet near $\varGamma \to 0$ for small times after $t_c$, $t>t_c$ where $\tau \equiv t-t_c\ll 1$. This is first explored by analytically continuing the Fourier series solution for the sheet shape $z$ obtained in (2.32) and (2.33) outside its domain of uniform convergence. Indeed, for $t< t_c$, the polylogarithm function $\textrm {Li}_{{5}/{2}}(w)$ contained in the summed form (2.33) with varying $\varGamma \in [0,2{\rm \pi} ]$ is evaluated inside the complex unit disk $\{w\in \mathbb {C},\,|w|<1\}$, where $\textrm {Li}_{{5}/{2}}({\cdot })$ remains a holomorphic function of $\varGamma$. When $t=t_c$, $\textrm {Li}_{{5}/{2}}({\cdot })$ is accessed along the unit circle, and its branch point $w=1$ is approached from both sides as $\varGamma \to 0$ in (2.33). As a result, the critical sheet shape $z_c$ loses analyticity at $\varGamma =0$, but remains continuous. For $t > t_c$, i.e. $\tau >0$, the asymptotic Fourier coefficients $\mathcal {A}_{n,0}$ given in (2.28) diverge as $n\to \infty$, and therefore the series solution breaks down. The summed form (2.33), originally obtained for $t < t_c$ can nonetheless be evaluated when $t>t_c$ as an analytic continuation because $\textrm {Li}_{{5}/{2}}(w)$ is in fact holomorphic for all $w\in \mathbb {C}\backslash (1,\infty )$. With $\lambda =1$ and $t= t_c+\tau$, this gives
where (2.31) is used, and
At the origin $z_{AC}(\varGamma \to 0^+,0)=0$ but since $J>1$ when $\tau >0$ the branch cut discontinuity along $w\in (1,\infty )$ of $\textrm {Li}_{{5}/{2}}(w)$ is crossed as $\varGamma \to 0^+$ in (2.39), leading to a discontinuous jump of the vortex-sheet profile at the origin. Specifically, this discontinuity at $\tau =0^+$ is given by $z_{AC}(0^+,0)=0$ when $\tau =0$ and
when $\tau >0$. Immediately after the critical time $t_c$, (2.41) is expanded for $\tau \ll 1$ as
suggesting that, when $\tau > 0$, the vortex sheet, which passes through the coordinate origin at $\tau = 0$, separates into two distinct branches (sheet tearing) whose ends or tips separate as $\tau ^ {3/2}$. This is shown in figure 3, where the analytically continued sheet shape is shown near the origin for two values of $\tau$.
In the following we put $b = |b|\exp (\textrm {i}\,\theta )$, $A = |A|\exp (\textrm {i}\,\chi )$, where $|b|$, $|A|$, $\theta$ and $\chi$ are known functions of $\epsilon$. We will refer interchangeably to variable sets $(\varGamma,\tau )$ and $(\eta,\tau )$ where
For clarity, we will refer to three regions defined for now by
(i) $\textrm {I}\quad \textrm {Outer:} \quad \eta \gg 1,\quad 2{\rm \pi} >\varGamma >0,\quad \varGamma \gg \tau$;
(ii) $\textrm {II}\quad \textrm {Intermediate:}\quad \eta = {O}(1),\quad \varGamma = {O}(\tau )$;
(iii) $\textrm {III}\quad \textrm {Inner:} \quad 0\le \eta <1,\quad 0\le \varGamma <\tau$.
Region I may be considered as one full spatial period of the infinite vortex sheet. The Moore approximation outlined above gives the region-I solution at $\tau =0^-$. Region II can be taken here as described by the initial condition (2.35) but more accurately as the region of validity of the region-II solution to be developed in § 5. There, it will be shown that analytic continuation of the region-II solution also exhibits sheet tearing. In Appendix B it is shown that sheet tearing noted above for the Moore solution and that for the region-II solution are somewhat different but agree in the limit $\epsilon \to 0$. Region III is to be discussed in § 6.
3. Intermediate problem: region II
A small-time solution in region II is sought as an expansion in a Taylor series about $\tau = 0$. This will be expressed by successive differentiation with respect to $\tau$ of the BR equation, leading to a formal expression for $\partial z^n(\varGamma,\tau \to 0)/\partial \tau ^n$ in terms of the partial Bell polynomials $B_{n,k}$. Evidence for the equivalence of the leading-order solution for the intermediate problem in region II and the small-time motion of the full initial-value vortex-sheet evolution in region I will be established by demonstrating close numerical agreement, when $\varGamma \to 0$, of the first two time derivatives for regions I and II for several $\epsilon$. These correspond to the initial angular velocity and angular acceleration, respectively, of the sheet motion near $\varGamma =0$. These terms will be seen to be singular as $\varGamma ^{-1/2}$ and $\varGamma ^{-3/2}$ respectively.
3.1. Expansion as a Taylor series
Since the integration domain is fixed and the integral is well defined in the Cauchy sense, the time derivative and the integral commute. We can then form the $n$th time derivative of (2.4) evaluated at $\tau =0$ as
The time derivative inside the integral can be expressed by making use of the Faà di Bruno (Johnson Reference Johnson2002) formula for the time derivative of the composite function $f(Z(t))$
where $f,Z$ are functions for which all derivatives are defined, the superscript written as $(k)$ denotes the $k$th time derivative and
are the partial Bell polynomials (Bell Reference Bell1927).
We apply (3.3) to evaluate the $n$th time derivative inside the integral in (3.1) by identifying $f=1/Z(t)$, where $Z(\tau ) \to z(\varGamma,\tau )\pm z(\hat {\varGamma },\tau )$ separately for each of the two terms, followed by summation. Use of (3.2) in (3.1) then gives an explicit expression
where $(\cdots \cdots +\cdots )$ denotes a second term with $z(\varGamma,0)-z(\hat {\varGamma },0)$ and associated time derivatives replaced by $z(\varGamma,0)+z(\hat {\varGamma },0)$, as per the second term in (2.4).
The first four time derivatives can be written as
where all are understood to be evaluated at $\tau =0$ with $z\to z(\varGamma,0)$, $\hat {z}\to z(\hat {\varGamma },0)$ and where the notation $[\cdots + \cdots ]$ indicates an addition of the preceding term with the minus sign replaced by plus. In principle, all derivatives can be evaluated sequentially. The first two terms of (2.35) are used in the kernel on the right-hand side of (3.5) and the integral evaluated to give $z^{(1)}(\varGamma,0)$. This can then be used in (3.6) to evaluate $z^{(2)}(\varGamma,0)$ and the process repeated. Each successive stage leads to increased complexity and, additionally, we cannot prove that all integrals converge.
3.2. Sheet velocity at $\tau =0$
Substituting (2.35) into the right-hand side of (3.5) and using the change of variables $\hat {\varGamma }=x\varGamma$ and ${\zeta } = A\sqrt {\varGamma }$ gives the complex velocity $u-\textrm {i}\,v$ at $\tau =0$ as
Expanding the integrand for small $|\zeta |$ gives
Using
then gives
The sheet velocity vanishes at the origin as $\varGamma \to 0$ at the critical time, but its angular velocity of order $z^{(1)}/\varGamma$ diverges as $O(1/\sqrt {\varGamma })$. So, for sufficiently small $0<\varGamma \ll 1$ and $\tau$, the vortex sheet admits the form
3.3. Sheet acceleration at $\tau =0$
Equation (3.6) can be written in full as
Using (3.13) at $\tau =0$ in (3.14) gives
Again expanding the integrand for $|\zeta | \ll 1$ and evaluating the CPV integrals gives, after some algebra,
Retaining only the leading-order term in (3.16), we can then extend (3.13) as
When the ordered expansions in $\tau \ll 1$, $\varGamma \ll 1$ are taken in (2.39), and then the $\epsilon \to 0$ asymptotes of both the resulting expression and (3.17) obtained, then asymptotic matching at least up to $O(\tau ^2)$ is found. A numerical example is discussed in Appendix B.
4. Numerical comparison at $\tau =0$
4.1. Periodic parametrization
The accuracy of (3.12) and (3.16) is tested by high-precision numerical evaluation of the BR integral (2.2). We periodically parameterize the interface shape $z(a)=x(a)+\textrm {i}y(a)$ and circulation $\varGamma (a)$ using ‘Lagrangian markers’, $a\in [0,2{\rm \pi} ]$, such that $x(a+2{\rm \pi} )=x(a)+2{\rm \pi}$, $y(a+2{\rm \pi} )=y(a)$ and $\varGamma (a+2{\rm \pi} )=\varGamma (a)+2{\rm \pi}$. The integral (2.2) can then be reduced to the finite domain
where $\sigma (a)=\textrm {d}\varGamma /\textrm {d}a$. Similarly, the acceleration integral becomes
where $z = z(a,t)$, $\hat {z}=z(\hat {a},t)$, $z^{(1)} = \partial z/\partial \tau$, $\hat {z}^{(1)} = \partial \hat {z}/\partial \tau$ and $\hat {\sigma }=\sigma (\hat {a})$. For numerical purposes, the introduction of interface markers $a$ allows discretized $\varGamma$ values to be distributed using a nonlinear monotonic stretching function of $a$, such as
4.2. Complex velocity near $\varGamma =0$
To evaluate the vortex-sheet velocity at the critical time $t = t_c (\tau =0)$, $z=z_c$ given by (2.34) is substituted into (4.1). With the choice of $\beta =0$, the principal-value integral can be evaluated directly for all $0\leq \varGamma \leq 2{\rm \pi}$ to arbitrary precision in a symbolic environment powered by mathematica®. A comparison between numerical evaluation of (4.1), denoted by $z^{(1)}_{num}$ and the intermediate-problem asymptotic approximation given by (3.12), and denoted presently by $z^{(1)}_{asy}$, is shown in figure 4(a), using $\epsilon =0.001, 0.01, 0.1$, where the $O(\sqrt {\varGamma })$ convergence of $z^{(1)}$ is clearly seen. Accuracy of the analytical estimate can be established for small $\varGamma$ in figure 4(b,c), where the relative error against $z^{(1)}_{num}$ decreases for both velocity components as $\varGamma \to 0$. For a given $\varGamma$, $z^{(1)}_{asy}$ associated with smaller initial perturbation $\epsilon$ has lower relative error dominated by the imaginary part.
4.3. Acceleration
The acceleration requires more care. An efficient numerical scheme to compute (4.2) was developed by first removing the singularity at $a= \hat {a}$ using the identity
Equation (4.2) can then be written as
where the integrand is continuous at $\hat {a}=a$ and takes the limit
Next, a uniform grid $\{a_k=2k{\rm \pi} /N\mid k=0,1,\ldots,N-1\}$ is formed, over which $z$ and $z^{(1)}$ at the critical time can be calculated with high precision using (2.34) and (4.1), respectively, as discussed in §4.2. This leads to the discrete sets $\{z_k\}$ and $\{z^{(1)}_k\}$, from which the $a$-derivatives of $z$ and $z^{(1)}$ are obtained using a Fourier method. Finally, (4.5) is summed with spectral convergence using the trapezoidal rule.
The difficulty with integrating (4.5) for small $\varGamma$ is that the desingularized kernel still exhibits sharp boundary layers near $\hat {a}=a$ and $\hat {a}=2{\rm \pi} -a$ when $a\ll 1$, which increase the truncation error. This is not surprising since $z^{(2)}$ is expected to be asymptotically divergent when $\varGamma \to 0$. The desingularized kernel, however, shows increasingly smooth behaviour for points $a_k$ further away from the origin. Therefore, the truncation error is mollified by choosing aggressive stretching $\beta =1$ in (4.3) and a large grid size $N=10^5$. The minimum $\varGamma$ obtained at $a=a_1$ is found at order $O(10^{-14})$, ensuring convergence of results obtained for the range of interest, e.g. $\sigma \gtrsim O(10^{-9})$. The comparison made between the analytic acceleration given by (3.16) and its numerical counterpart in figure 4(d–f) shows satisfactory agreement between the two for small $\varGamma$. Relative errors comparable to those seen for the asymptotic velocities are also observed here, where the accuracy of the analytic estimate suffers significantly when $\epsilon$ and $\varGamma$ are increased.
5. Further Taylor-series expansion in region II
We now extend the analysis of § 3.1 by considering a series solution of the form
Use of (5.1) in (3.4) gives a triple series resulting in, at order $n$, a combinatorial explosion which is presently intractable. In order to make progress, simplifications are implemented. The first is truncation of the inner series in (5.1) at $m=1$ to give a series that can be rearranged and expressed as
with $C_0,C_1,C_2$ known and $C_n,n>2$ to be determined. The first four terms of (5.2) agree with (3.17), the first two of which are the region-II initial condition, while the second two were shown in § 4 to give good agreement with numerical calculations. We now proceed to extend the preliminary solution (3.17) using (5.2). This will later be shown to provide a good approximate solution of an appropriate linearized BR equation in region II.
When (5.2) is used in (3.4), at each order $\tau ^n$, the Bell-polynomial expansion produces a series in $\varGamma$ in a sum over negative powers for $n\ge 2$: ${\varGamma ^{3/2-n},\, \varGamma ^{2-n},\varGamma ^{5/2-n},\ldots}$ At order $n$ we consider only the most singular term of order $\varGamma ^{3/2-n}$, $\varGamma \to 0$, which is expected to be dominant in the subsequent analysis. Further, an induction-based analysis of the Bell-polynomial expansion (3.4) indicates that, at the order of the $n$th $\tau$-derivative, the most singular term of order $\varGamma ^{3/2-n}$ is contained in the $k=2$ term displayed as the right-most term with denominator $(z(\gamma,0)-z({\hat {\varGamma }}))^2$ in the four $\tau$-derivatives given in (3.5)–(3.8). For general $n\ge 2$, this can be expressed as
where, for clarity and until otherwise needed, the ‘II’ subscript has been suppressed. Equation (5.3) is the basis of further analysis.
5.1. Recurrence relations
From (5.2) we can write
In the following, we make the further approximation that, in the denominator of (5.3), we can write $z(\varGamma,0)- z(\hat {\varGamma },0) \approx b(\varGamma -\hat {\varGamma })$. This can be justified by observing that the additional $\varGamma ^{3/2}$ term in the initial condition (2.35) will, for small $\varGamma$, produce a first correction consisting of an additional power $\varGamma ^{1/2}$, as was seen in the earlier calculation of the initial velocity and acceleration terms. At the time derivative of $o(n)$, this will not affect the most singular contribution proportional to $\varGamma ^{3/2-n}$, but will contribute to less singular powers of $\varGamma$.
Using (5.4) on the left-hand side of (5.3) and its complex conjugate with $n\to n-1$ in the integrand on the right-hand side, leads to cancellation on both sides of factors $\varGamma ^{3/2-n}$ and, after some algebra, to the recurrence relation
The integral exists as an CPV integral for ${\tfrac {1}{2}} < n < {\tfrac {9}{2}}$. Presently, we consider only the finite part, taken as its analytic continuation for $n>{\tfrac {9}{2}}$. For real integer $n$ this is
Now make the change of variables
so that (5.2) can be written as
The similarity variable $\eta$ will be used throughout the sequel. Using (5.5) and $n\to n+1$, the recurrence relation for $K_n$ is
which, with iteration, can be written as
where we note that the recursion coefficient is real.
The above suggests splitting (5.8) into two series as
with $\hat {K}_0=1$, $\hat {K}_1=3/4$ and the $\hat {K}_n$ satisfy the same recurrence relation as do the $K_n$, namely (5.10). Solution of (5.10) gives for $n$ even and odd, respectively,
where ${\hat {\varGamma }}$ here denotes the complete gamma function, to be distinguished from the circulation variable $\varGamma$. Substitution of (5.13a,b) into (5.12a,b) shows that both series converge for $\eta >1/2$ or for $\varGamma > \tau /(2\,|b|^2)$. Hence, $\eta =1/2$ is now defined as the inner boundary of region II. For finite $\tau >0$ (5.11) cannot be extended to $\varGamma \to 0$.
5.2. Closed-form solution
The two series (5.12a,b) with $\hat {K}_n$ given by (5.13a,b) can be summed analytically, and when multiplied by $\eta ^{3/2}$ as suggested by (5.11), can be expressed as
The functions $Q_1(\eta ),\,Q_2(\eta )$ have fractional-power branch-point singularities in the complex $\eta$-plane at $\eta = \pm \textrm {i}/2$ and at $\eta =0$, confirming the radius of convergence of the series (5.12a,b). When used in (5.11), this gives
In the complex $\varGamma$ plane, with $\tau$ considered a real, positive parameter, the right-hand side of (5.16) has fractional-power branch points at $\varGamma = \pm \textrm {i}\tau /(2\,|b|^2)$ and at $\varGamma \to \infty$. When $\tau$ increases from $\tau =0^+$, the branch point at $\varGamma =0$ at $\tau =0$ splits into two points, each moving away from $\varGamma =0$ along the positive and negative imaginary axes, respectively, and $z(\varGamma,\tau )$ given by (5.16) is analytic at $\varGamma =0$.
5.3. Analytic continuation
We take (5.16) as the region-II solution for $\eta \ge {\tfrac {1}{2}}$, which corresponds to $\varGamma \ge \tau /(2\,|b|^2)$. This equation can nonetheless can be evaluated inside the circle $|\eta | < {\tfrac {1}{2}}$ in the complex $\eta$-plane or $|\varGamma |< \tau /(2|b|^2)$, and therefore gives the analytic continuation of (5.12a,b)–(5.13a,b) in this region. This function is finite valued and analytic on the real line in $0 \le \varGamma <\infty$. Its properties and relation to (2.39) when $\varGamma \to 0$ at fixed $\tau$ are discussed Appendix B where it is shown that this solution gives sheet rupture and separation into two branches at $\varGamma = 0^+$ for $\tau >0$.
In $(\eta,\tau )$ variables we write (5.16) as
where, for convenience, the $\exp ({\textrm {i}\,\theta })$ term is included in the definition of $\varOmega _{II}$ to denote rotation to the sheet tangent at $z=0$, allowing for clarity in the discussion of the $w_0(\eta )$ and $w_1(\eta )$ functions. Expressing the first few terms of the series ((5.12a,b)–(5.13a,b)) gives in $(\varGamma -\tau )$ variables
the first four terms of which agree with (3.17). In terms of $(\eta -\tau )$ variables, the first few terms of $w_1(\eta )$ are
6. Inner solution: region III
6.1. Separated algebraic-spiral structure
Here, we discuss an inner region III defined presently by $0 \le \varGamma < O( \tau /(2\,|b|^2))$, which exists only for $\tau >0$. The analytic behaviour of (5.16) when $\eta \to 0$ at fixed $\tau$, equivalent to $\varGamma \to 0$, is discussed in Appendix B, where it is shown that, provided $A_r \ne A_i$ (recall $A= A_r + \textrm {i}\,A_i$), this shows sheet end separation at $\tau = 0^+$ into two distinct branches whose endpoints move away from the origin as $\tau ^{3/2}$. This is non-physical because a vortex sheet with a free-end point is expected to roll up. Further, in Appendix C, it is shown that (5.18a,b) is an accurate approximate solution of a linearized BR equation for $\eta >1/2$ but generally fails when $\eta < 1/2$. An exceptional case is an initial condition (2.35) with $A_r = A_i$, where (5.16) then joins the origin smoothly. An alternative structure is required for the general case $A_r \ne A_i$, which is now discussed.
Equation (5.17a,b) suggests a general ansatz for an inner solution of the form
Truncating (6.2) at $\tau ^{1/2}$ inside the outer brackets and substituting into (2.4) gives
where $\omega _0 = \omega _0(\eta ),\omega _1 =\omega _1(\eta ),\omega _0' = \omega _0(\eta '),\omega _1' = \omega _1(\eta ')$. Equation (6.3) is a single equation for two functions $\omega _0(\eta ),\omega _1(\eta )$. For $\tau \ll 1$ (6.2) is considered as the leading order in an expansion of an inner solution in the small parameter $\tau ^{{1}/{2}}/|b|$. Further terms $\tau /|b|^2\,\omega _2(\eta ), \tau ^{{3}/{2}}/|b|^3 \omega _3(\eta )\,\ldots$ can be added but are not considered presently.
Expanding the integrand of the right-hand side of (6.3) and equating powers in $\tau ^{{1}/{2}}/|b|$ gives, to order one and $\tau ^{{1}/{2}}/|b|$, respectively,
Equations (6.4) and (6.5) are nonlinear and linear singular integral-differential equations for $\omega _0(\eta )$, $\omega _1(\eta )$, respectively. Equation (6.4) is a special case with $p=1$ of a more general similarity integro-differential equation with dimensionless parameter $p$ such that $\omega _0(\eta )\sim \eta ^{1/p}$, $\eta \to \infty$ (Pullin & Phillips Reference Pullin and Phillips1981; Pullin Reference Pullin1989). Boundary conditions with $p=1$ for $\eta \to \infty$ are given by matching to the intermediate solution in region II giving $\omega _0\sim \eta$ (first of (5.18a,b)), $\omega _1 \sim A \eta ^{{3}/{2}}$ (5.20). Equation (6.4) has the exact solution $\omega _0= \eta$ matching (5.18a,b). This corresponds to a uniform, flat vortex sheet. The status of (5.18a,b) as an approximate solution to (6.5) in region II but with $(\omega _0,\omega _1)$ replaced by $(w_0, w_1)$, respectively, is discussed in Appendix C. It is concluded that, with $\omega _0= \eta$, smooth solutions of (6.5) that match the region-II solution and that are bounded when $\eta \to 0$, do not exist. Hence, $\omega _0=\eta$ cannot provide a basis for an admissible region-III solution.
6.2. Numerical solution: zeroth order
Pullin (Reference Pullin1989) reported numerical evidence that an alternative solution to (6.4) exists, satisfying $\omega _0\sim \eta$, and taking the form of two distinct and separated vortex cores, each with centres of counter-clockwise roll-up removed from the origin. We denote these roll-up points in the $\eta$ plane as $\omega _{0,c} = \omega _0(\eta =0)$ and $-\omega _0(\eta =0)$ for the right-hand ($\mathrm {Re}(\omega _{0,c})>0$) and left-hand ($\mathrm {Re}(\omega _{0,c})<0$) sheet sections, respectively, consistent with the required symmetry of the overall solution. Evidence was also provided (Pullin & Phillips Reference Pullin and Phillips1981; Pullin Reference Pullin1989) that solutions with $\omega _{0,c}=0$, and satisfying the same far-field boundary condition, which corresponds to spiral roll up about the origin, do not exist.
Non-existence was inferred in the sense that, with fixed $\omega _{0,c}=0$, while double-spiral solutions were found for $p<1$, in the limit $p\to 1^-$ from below their spatial scale in the $\omega$ plane shrunk to zero. In contrast, numerical solutions on a different branch showed smooth variation as $p$ passed through unity from above ($p \to 1^+$), giving a particular solution with $p=1$ consisting of separated algebraic spirals of finite scale and with $\omega _{0,c}\ne 0$. This vortex structure is consistent with end-sheet separation of the $\eta \to 0$ continuation of the region-II solution, while exhibiting free-end roll up expected on physical grounds. The above does not rule out that a double-spiral solution to (6.4) may exist either in isolation or on an undiscovered solution branch, but, despite searching, none has been found. We will interpret the separated-sheet solution as the only available and admissible zeroth-order, region-III solution that allows matching to the region-II solution.
6.2.1. Numerical method
The calculations of Pullin (Reference Pullin1989) have been reproduced to give a solution for $\omega _0$. The original code and specific results are no longer available, and so independent coding was implemented using the same method. This is now summarized briefly. The real line $\eta \ge 0$, $(0,\infty )$ is divided into three contiguous sections $(0,\eta _N),(\eta _N,\eta _0)$ and $(\eta _0,\infty )$. The vortex sheet in these sub-domains is represented by three parts: (i) an inner section in $0\le \eta \le \eta _N$ modelled by a point vortex of strength $\eta _N$ located at $\omega _{0,c}=\omega _{0,N+1}$, (ii) an intermediate, continuous section in $\eta _N\le \eta \le \eta _0$ represented by $N$ straight segments, $(\omega _{0,j}-\omega _{0,j-1}), j = 1,\ldots,N$ and (iii) an outer section $\eta _0\le \eta \le \infty$ modelled by $\omega _0(\eta ) = \eta$. In part (i), $\omega _{0,N}$ is joined to $\omega _{0,N+1}$ by a cut in the $\omega$-plane in order to define a single-valued complex velocity potential. The unknowns are the $2N+2$ quantities given by the real and imaginary parts of $\omega _{0,j}$, $j=1\ldots N+1$.
In (ii), a finite-difference form of (6.4) is satisfied at the midpoint of each segment using a two-point rule for derivatives on the left-hand side, a trapezoidal rule for the Cauchy principal-value integral contributed from (ii), and a point-vortex contribution from (i). The contribution from (iii) can be evaluated analytically. This gives $2 N$ nonlinear equations. The trapezoidal rule incurs error owing to neglected logarithmic corrections arising from proximity of a point on the sheet to nearby segments on adjacent spiral turns. This is somewhat relieved by cancellation since, as will be seen, most points lie inside a tightly wound spiral with many individual turns on either side. In (i), pointwise representation of (6.4) is replaced by an integrated approximation over $0 \le \eta \le \eta _N$, giving two additional equations. The $2(N+1)$ equations are solved by a Newton–Raphson method with analytical calculation of the Jacobian.
Numerical solutions produce a tightly wound, algebraic spiral with many turns. Related BR solutions driven by singularity formation have been reported by Pullin & Phillips (Reference Pullin and Phillips1981), Pullin (Reference Pullin1989) and Pullin & Sader (Reference Pullin and Sader2021), using a similar numerical method. Convergence with respect to $N$ and the number of spiral turns captured was documented by Pullin & Sader (Reference Pullin and Sader2021). Present solutions use $N=1016$ to $N= 2816$ with $\eta _0 = 12.5$ chosen to be remote from where the solution shows significant deviation from $\omega _0 = \eta$. Testing with different values of $\eta _0$ showed negligible difference provided $\eta _0 \ge 2$. With $\eta _0=12.5$ numerical solutions near this value agree with the asymptotic form $\omega _0=\eta$ to order $10^{-6}$ for $\mathrm {Re}[\omega _0]$ and order $10^{-4}$ for $\mathrm {Im}[\omega _0]$. Some parameters for numerical solutions are listed in table 1. The portrait of the parametric solution in the $\omega _0$ plane is shown in figure 5 for $N = 2816$. Solutions with $N = 1916-2216$ are indistinguishable to plotting accuracy. The shape closely resembles the same flow as in Pullin (Reference Pullin1989), who used $N= 170$ with a few turns. In the present scaling his solutions show $\omega _{0,N+1} \approx (0.22+ 0.075\,\textrm {i})$ compared with the more accurate values listed in table 1. The non-zero value $\omega _{0,N+1}$ for the right-hand sheet branch can be interpreted as a balance between the Biot–Savart induction of the far-field circulation distribution and that of the compact, left-hand vortex core.
6.2.2. Inner structure of zeroth-order solution
An image of the inner portion of the sheet geometry of the right-hand core is shown by the solid line in figure 6, which indicates a tightly wound, almost circular spiral with structure strongly dominated by the Biot–Savart induction at a point $\omega _0(\eta )$, provided by the cumulative sheet circulation ($\eta$) that lies inside a circle that passes through $\omega _0(\eta )$ and is centred on $\omega _{0,c}=\omega _{0,N+1}$. This model approximation goes back to Kaden (Reference Kaden1931). It does not include the induction effects of the remainder of the spatial circulation distribution that lies outside this circle and that includes the left-hand rolled-up vortex. The inner portion can then be analysed by replacing (6.4) by the model equation
where $\tilde {\omega }_0 = \omega _0 - \omega _{0,N+1}$, with exact solution for $\tilde {\omega }_0(\eta )$ and corresponding zeroth-order $\tilde {z}_{III}(\varGamma,\tau )$ given by
where $\alpha,\,\phi _0$ are arbitrary, real parameters. Sohn (Reference Sohn2016) shows that (6.7) satisfies (6.4) to a good approximation. Here, we patch (6.7) to an arbitrary point on the inner part of our numerical solution for $\omega _0$. We use $\eta = 0.02389, \omega _0= 0.2197-0.07855\textrm {i}$, which is well inside the rolled-up sheet. Other choices show numerically similar results. This gives two equations for $\alpha,\phi _0$ with numerical solution $\alpha =0.2397\,,\phi _0 = 0.7886$. In figure 6, the patch point is marked by a black dot. The numerical solution with $N=2816$ contains $15$ sheet turns inside the black dot. Equation (6.7) is plotted as a dashed line in the range $0.008\le \eta \le 0.023886$. The agreement with the numerical solution is satisfactory with systematic differences attributable to a slight ellipticity in the numerical solution produced by the combined straining effect of the vorticity structure that lies outside the inner portion, and that includes the other vortex core.
The inner structure of the vortex core in the $\omega$ plane can be estimated using (6.7) where the radius of the spiral decreases as $|\omega _0|\sim \alpha \eta$ while its angle relative to an arbitrary datum increases as $\theta _s \sim 1/(2\,{\rm \pi} \alpha ^2 \eta )$. Hence, the ratio of local spiral turn spacing to local radius decreases as $\delta |\widetilde {\omega _0}|/|\widetilde {\omega _0}| \sim \theta _s^{-1}$, which decreases rapidly as $\theta _s$ increases. This variation, a result of the essential-singularity (in a complex $\eta$ plane) character of the $\eta \to 0$ limit solution, suggests potential difficulties in time-domain numerical solutions of the BR initial-value problem from $\tau = 0^- \to \tau =0^+$. In the $z$-plane, (6.8) shows that the smoothed azimuthal velocity distribution is $u_\theta =1/(2\,{\rm \pi} \alpha |b|)$ (with jumps across each sheet turn) and is independent of radius and $\tau$ while the smoothed vorticity distribution is singular as $u_\theta /r$, where $r$ is the radius to the roll-up centre. If one considers a finite portion of the right side ($x>0$), almost flat vortex sheet bounded at $\tau =0$ by the origin (singular point) and an arclength $x=S$ along the sheet, then this will contain circulation $\varGamma (S) \approx S/|b|$. This portion will roll up to constant radius $r(S) =\alpha \,S$. The quantity $\beta = 1/(2\alpha )$ is known as the Betz constant (see Moore & Saffman Reference Moore and Saffman1973) with value $\beta \approx 2.08$.
6.3. Numerical solution: linearized first-order equation
The numerical solution for $\omega _0(\eta )$ can now be substituted into (6.5) to give an equation for $\omega _1(\eta )$. This was solved presently with essentially the same numerical method described above for $\omega _0$, except that Newton iteration is not required owing to linearity. At $\eta =\eta _0$, the boundary condition is $\omega _1(\eta ) = w_1(\eta )$ given by the second of (5.18a,b), allowing effective matching of the region-III solution with that in region II. The $\eta \to 0$ limit for each solution is represented by the isolated point $\omega _{1,N+1}$ which is determined as part of the solution. Because the constant $A$ is a function of $\lambda _1\epsilon$, then solutions to the linearized equation will be a one-parameter family. An asymptotic limit solution will be discussed subsequently. Testing of the numerical method for the solution of (6.5) is described in Appendix C where it is shown that only first-order convergence is achieved. With $N$ large, this is considered satisfactory.
Solution portraits in the $\omega _1$ plane are shown in figure 7(a–c). These panels show self-intersection, which would be unacceptable for the full-sheet shape but not for the linearized solution in the form shown. Figure 8 gives the variation of both $\mathrm {Re}[\omega _1(\eta )]$ and $\mathrm {Im}[\omega _1(\eta )]$ for several values of $\epsilon$ in the range $10^{-2}-10^{-10}$, with $\lambda _1=1$. Dashed lines are the region-II solutions given by the second of (5.18a,b). Numerical solutions show oscillation about $\omega _{1,N+1}$. In fact the plots in figure 8 suggest that, as $\epsilon$ is reduced, $\omega _{1,N+1}$ appears to approach close to the limit of the region-II solution for $w_1(\eta \to 0)$. In turn this suggests the existence of an $\epsilon \to 0$, asymptotic solution.
This can be explored by taking the limit $\epsilon \to 0$ of $(b(\epsilon ),A(\epsilon ))$ defined by (2.36a,b), giving
The $\epsilon \to 0$ asymptotic form of (5.18a,b) then becomes
This gives the ${\hat {\varepsilon }}\to 0$ limit solution $\omega _1(\eta ) = 2{\hat {\varepsilon }}\hat {\omega }_1(\eta )/3$ which scales $\epsilon$ out of the linear problem. Equation (6.5), (with $\omega _0(\eta )$ given by the numerical zeroth-order solution) was solved numerically for $\hat {\omega }_1(\eta )$ with boundary condition given by $\hat {\omega }_1(\eta _0) = \hat {w}_1(\eta =\eta _0)$. The portrait of the solution in the $\hat {\omega }_1$ plane, obtained with $N= 2816$, $\eta _0= 12.5$, is shown in figure 7(d). The calculated value of the end or centre point is $\hat {\omega }_{(1,N+1)} = -0.5607+ 0.5278\textrm {i}$ compared with $\hat {w}_1(\eta \to 0) = -0.5+ 0.5\textrm {i}$, giving limit ratios $\mathrm {Re}[\hat {\omega }_{(1,N+1)}]/\mathrm {Re}[ \hat {w}_1(0)] = 1.121$ and $\mathrm {Im}[\hat {\omega }_{(1,N+1)}]/\mathrm {Im}[ \hat {w}_1(0)] = 1.056$. Figure 9 shows the ratios $\mathrm {Re}[\omega _{(1,N+1)}]/\mathrm {Re}[ w_1(0,\epsilon )]$ and $\mathrm {Im}[\omega _{(1,N+1)}]/\mathrm {Im}[ w_1(0,\epsilon )]$ obtained from numerical solutions for finite $\epsilon$. The latter shows rapid convergence towards the asymptotic limit, while the convergence of the ratio of the real components is much slower.
6.4. Composite solutions
We illustrate the composite solutions for $\tau >0$ with $\epsilon = 10^{-3},\lambda _1=1$. These are calculated using (6.1) with (6.2) truncated at order $\tau ^{1/2}$. It is noted that the basic solutions $\omega _0,\omega _1$ described earlier were obtained, for convenience, in a reference frame where the zeroth order is $w_0=\eta$. In order to properly match the region-II solution given in (5.18a,b), (6.2) shows that these must be rotated anticlockwise by $\theta (\epsilon )= \arg (b(\epsilon ))$
Figure 10 shows two images of the sheet evolution at in the $\varOmega$ plane at $\tau = 0.001, 0.01$. The dashed line shows the zeroth-order sheet shape represented by $\omega _0$. The order $\tau ^{1/2}$ correction by the outer flow appears to shift the sheet profile, relative to its zeroth-order position, towards the origin. Since this will be active for both sheet branches, the overall effect is to convect the two sheets towards each other. Successively expanded views of the double-sheet structure are shown in figures 11 and 12 at $\tau = 0.01$, where the intermediate sheet profile $\varOmega _{II}(\eta )$, defined in (5.18a,b), is also shown as dashed lines. Figure 11 displays the scale of the separated spirals and also the matching of the inner region-III solution onto the region-II solution. Figure 12 illustrates the spatial scale of the rolled-up spiral in comparison with the large radius of curvature of the region-II solution. Time-domain images can be obtained by multiplying the $\varOmega$ solutions by $\tau /|b|$. These are shown in figure 13 for $\tau = 0.0005 ,0.001,0.0025 ,0.005,0.0075,0.01,0.0125$. The effect of the region-II correction is discernable but is small.
7. Conclusion
We have considered the evolution of a periodically disturbed vortex sheet during the nonlinear stage of spatially periodic, KH instability for small times $\tau = t-t_c$ subsequent to the formation of a curvature singularity as analysed by Moore (Reference Moore1979), for an initial sheet perturbation defined by the small parameter $\epsilon \ll 1$. The time-continued solution is itself singular, comprising three regions in the space of the variables $\eta \sim \varGamma /\tau$ and $\tau$. The outer region I denotes the full extent of a single spatial period of the vortex sheet whose shape and strength at $\tau =0$ are defined by an extension of Moore's asymptotic solution for a general initial condition, giving a closed expression for the sheet evolution up to $t=t_c$ in terms of the polylogarithm function. Its analytic continuation beyond the radius of convergence of Moore's Fourier series solution corresponds to $\tau >0$. This shows sheet rupture with end separation proportional to $\tau ^{3/2}$.
An intermediate domain (region II) is identified as a vortex-sheet evolution with initial condition at $\tau =0$ given by the first two terms of the expansion of the region-I solution near the singular point $z=0$ in the Lagrangian circulation variable $\varGamma$. This contains the essential elements of the singularity structure as a $\varGamma ^{3/2}$ term added to the otherwise locally flat sheet, including known quantitative dependence of coefficients on the small-amplitude parameter $\epsilon$ contained in the region-I initial condition. This is done for a general initial sheet-shape perturbation. For small $\tau >0$, this initial condition is evolved to $\tau >0^+$ using a Taylor-series expansion based on the dominant term arising from time-wise differentiation of the governing BR equation. The first two terms corresponding to the velocity and acceleration show good agreement with Biot–Savart numerical calculations using the full region-I solution. Both terms contain singular features.
The Taylor series has a finite inner radius of convergence in the complex $\eta$ plane, but can be summed to provide a closed-form expression which can be analytically continued to $\eta \to 0$, equivalent to continuation to zero circulation $\varGamma$ at fixed $\tau >0$. This also shows sheet tearing or jumping at $\tau =0^+$, in qualitative agreement with the behaviour of the Moore (Reference Moore1979) solution, analytically continued to $\tau >0$. The closed-form expression is shown in Appendix C to be an accurate approximate solution of the linearized BR equation in region II, but generally fails for $\eta \to 0$. With one exceptional case, a solution of the linearized BR equation that both matches the region-II solution, and also is bounded near the point of original singularity formation, could not be found. This exceptional case corresponds to $\arg (A)= {\rm \pi}/4,\, 5{\rm \pi} /4$ where $A(\epsilon )$ is the strength of the singular part of the sheet-shape function at the critical time. This case is not encountered within the class of initial conditions with $\mathrm {Re}[A]>0$.
The tearing scenarios found from both the analytically continued region-I and region-II solutions show flat sheet ends, without vortex-sheet roll up. While considered non-physical in itself, this points to the existence of a distinct, inner region III with a local solution structure that incorporates an analogue tearing event. An inner solution in a region III is explored with analytical form suggested by the structure of the region-II solution. In the $z$-plane, but expressed in terms of $(\eta,\tau )$ variables, this is proportional to $\tau$ times the first two terms of a perturbation series in $\tau ^{1/2}$ with coefficients defined by functions determined respectively and recursively by solutions to a nonlinear zeroth-order ($O(\tau ^{0})$) and a first-order ($O(\tau ^{1/2})$) linearized BR equation, both in similarity form with independent variable $\eta$. The region-III solution exists only for $\tau >0$.
Numerical solutions of both region-III similarity equations are obtained that effectively match the corresponding region-II solutions at $\eta \to \infty$. The zeroth-order solution is independent of $\epsilon$, showing an antisymmetric pair of separated, algebraic-spiral vortex sheets of well-known form, each with its own roll-up centre. Tearing or rupture at $\tau =0^+$ is interpreted as the response of the vortex sheet, which is continuous at the origin at $\tau =0^-$, to the onset of singularity formation which prevents its smooth and contiguous continuation for $\tau >0$. It is argued that an inner, zeroth-order solution based on finite-sized self-similar, double spiral roll-up about a common centre at the origin does not exist, although no proof is given. At first order $\tau ^{1/2}$, numerical solutions of the linearized BR equation exhibit oscillatory behaviour forced by the spiral form of the zeroth-order solution. These show dependence on $\epsilon$. An asymptotic form is obtained with explicit analytic dependence on $\epsilon$.
In summary, at $\tau =0^-$ the vortex sheet is continuous but with singular curvature (Moore Reference Moore1979), angular velocity, acceleration and other higher-order derivatives at the local centre of antisymmetry. At each point of infinite curvature at $t=t_c$ ($\tau =0$), this singular behaviour generates local sheet rupture at $\tau =0^+$ into two distinct branches at the central singular point with separation distance of the branch endpoints, or roll-up centres, scaling as $\tau$ together with a non-local, $\tau ^{3/2}$ correction. In the complex $\varGamma$ plane this can be interpreted as the spontaneous appearance of a pair of isolated, essential singularities. In the extended complex $z$-plane, with a periodic initial perturbation, one singular point exists in each spatial period. The vortex sheet, continuous and of infinite streamwise extent for $t\le t_c$, fractures at $\tau =0^+$ into an infinite array of identical vortex sheets, each with two algebraic-spiral, roll-up endpoints generated at neighbouring singular points. While no attempt has been made presently to determine the character of the complex velocity for $t>t_c$, it is plausible, at least for sufficiently small $\tau$, that this is analytic on a Riemann surface defined on a strip of width $2{\rm \pi}$ streamwise and of infinite extent in the $\pm \mathrm {Im}[z]$ direction. The vortex sheet connecting the two roll-up centres would define the locus of a branch cut. Without further analysis it is not clear if the spiral centres themselves are standard branch points or conform to a different singular form.
Finally, we remark that the present solution provides a classical vortex-sheet evolution from a singular initial condition at $t=t_c$. Its uniqueness has not been demonstrated. A further interesting question is whether this solution is, in some properly defined sense, the unique limit solution of either the Euler or Navier–Stokes equations at some fixed $t>t_c$ when an appropriate regularizing parameter tends to zero. This remains open.
Acknowledgement
The authors thank G.R. Baker for helpful discussions.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Uniformly convergent Fourier series
In order to extend the Moore (Reference Moore1979) analysis to a general initial condition, we briefly derive (2.28). Following Moore (Reference Moore1979), we first construct a generating function $G(x,t)\in \mathbb {C}$, with real $x$ and $t$, as follows:
where $\mathcal {A}_{n,0}$ are the leading-order Fourier coefficients defined in (2.12). As a result, the ODE system (2.14) for $\mathcal {A}_{n,0}$ is equivalent to
where the overline denotes complex conjugate. To obtain the asymptotic form of $G(x,t)$ for large $t$, (2.18) is substituted into (A1), giving
where
and $\lambda _{j_r}$ are the constants found in (2.18). However, as discussed in § 2.2.1, expansion (A3) leads to non-uniformly valid coefficients $\mathcal {A}_{n,0}$.
To resolve this issue, a strained time $s(t)$ is introduced
where constant coefficients $\alpha _n$ are to be determined. Integrating (A5) leads to
Consequently, (A3) can be rearranged as
where
and $h_r$ can be obtained by matching (A7) with (A3). As $\mu \to 0$, one has for example,
where $\lambda _{1,0}$, $\lambda _{2,0}$ and $\lambda _{2,1}$ are given in (2.17).
Alternatively, directly substituting (A7) into (A2) gives a recursive ODE system for $h_r$, where the first two equations read,
Solutions of (A10) subject to initial values (A9a,b) are
where $g(\mu )$ is given by (2.23), and
The constant $\alpha _1$ is determined by expanding (A12) around the branch point $\mu =-1/(\lambda _{1} e)$ and requiring that the most singular term in the expansion is no worse than that of (A11), i.e. $O(-1/e+\lambda _{1}\mu )^{3/2}$. This produces
which is non-zero in general. For the Moore case of $\phi ={\rm \pi} /2$, we recover $\alpha _1 = 2$. Finally, substituting (A11) and (A12) into (A1) leads to the desired result (2.28).
Appendix B. Analytic continuation of region-II solution within region III
B.1. Inner expansion of analytically continued region-II solution
The right-hand side of $w_1(\eta )$ given by the second equation in (5.18a,b) can be expanded about $\eta =0$ to give
which then gives, for the inner expansion of $z_{II}(\varGamma,\tau )$
When $\varGamma \to 0^+$ we find that
This result states that, when $\tau >0^+$, the vortex sheet, which passes through the coordinate origin at $\tau =0^-$, separates into two distinct branches whose ends or tips separate as $\tau ^{3/2}$. An exceptional case is $A_r = A_i$. A similar result was also obtained as (2.42) from the analytic continuation of the Moore solution. The real and imaginary parts of both expressions are plotted vs $\epsilon$ in figure 14. When $\epsilon \to 0$ both expressions can be shown to be asymptotic to
Agreement between the two expressions is better for the imaginary part. The sheet separation at $\tau = 0^+$ appears to be an inescapable consequence of the analytic continuation of Moore's solution (2.33) which we expect to be asymptotically accurate when $\epsilon \to 0$.
In fact the correspondence of two solutions can be taken further. If from (2.39) we form $z_{AC}(\varGamma,\tau )-z_{AC}(\varGamma,0)$, then this will give the change in sheet shape from its initial condition at time $\tau$. We can do the same by subtracting (2.35) from (5.16). Real and imaginary parts of these difference expressions are plotted in figure 15 for $\tau = 10^{-4}, \epsilon = 10^{-3}$, where it can be seen that relative differences are generally $O(1)$. In figure 16 using $\epsilon = 10^{-20}$, differences are smaller especially for the imaginary part. The reasons for this are not clear but are perhaps that first, that the assumption $\epsilon \ll 1$ is built into the Moore solution but not explicitly assumed in the Taylor-series expansion leading to (B2). Second, the main results of the Moore solution, surprisingly, follow from linear equations given by (2.14) while (5.16) is expected to be valid only for small $\tau$. We note that the series form of (5.16) is convergent only for $\varGamma > \tau /(2|b(\epsilon )|^2)$. With $\epsilon = 10^{-3}, 10^{-20}$ this gives $\varGamma > 7.46\times 10^{-5}$, $\varGamma > 5.25\times 10^{-5}$ respectively. These values both lie within the range of $\varGamma$ plotted in figures 15 and 16.
B.2. Sheet strength, $\varGamma \to 0$, $\tau >0$
From (B2) we find that to $O(\eta )$
The vortex-sheet strength $\gamma (\varGamma,\tau ) = |\partial z_{II}/\partial \varGamma |^{-1}$ for small $\eta$ can be estimated. At the right-hand sheet tip when $\eta \to 0$
This is finite for $\tau =0^+$ with an order $\tau ^{1/2}$ correction for $\tau >0^+$.
Appendix C. Linearized BR equation in region II
We investigate $w_1(\eta )$ given by the second of (5.18a,b) as a solution of (6.5), written here with $w_0 = \eta$ as
Using that
we obtain
The double bar indicates a generalized Cauchy principal-value integral interpreted as the average of values obtained when approaching $\eta$ on the real line from above and below, in the complex $\eta$ plane.
To proceed, we employ an approach utilized by Cowley et al. (Reference Cowley, Baker and Tanveer1999, § 5). First write $w_1 = w_r +\textrm {i}\,w_i$ and introduce the splitting
which, when used in (C3), gives real and imaginary parts respectively as
With the constant $A = A_r+\textrm {i}\,A_i$, (5.18a,b) and (C4a,b) gives test solutions as
which can be expressed as
With these expressions the integrals on the right-hand sides of both (C6) and (C7) converge at both $\eta ' =0,\infty$. Cowley et al. (Reference Cowley, Baker and Tanveer1999) used a Fourier-transform method to show that $F_1(\eta )$ is an exact solution of (C6). This has been verified presently both using the sine transform on $0\le \eta <\infty$ and also numerically. The function $F_2(\eta )$ is not an exact solution of (C7). However, a series solution of (C7) in inverse powers of $\eta$ can be constructed for large $\eta$. This is straightforward and details are omitted. The series is convergent for $\eta >{\tfrac {1}{2}}$ and sums to $F_2(\eta )$, which is then presently interpreted as an accurate approximate solution to (C7) for $\eta > {\tfrac {1}{2}}$ (this is verified numerically below) and therefore an asymptotic solution when $\eta \to \infty$. We conclude that $w_1(\eta )$ given by the second of (5.18a,b) is generally a very good approximate solution of the linearized BR equation (C3) for $\eta >{\tfrac {1}{2}}$ in region II, an asymptotic solution when $\eta \to \infty$ but fails for $\eta \to 0$. The exceptional case is $A_r = A_i$ when $g^-_s =0$ and $g^+_s$ gives an exact solution for all $\eta$.
The solution $F_1(\eta )$ for (C6) can be used as a test for the numerical method used for the solution of the linearized BR equation described in § 6.3. This was implemented as a numerical solution of (C6) in real variables on the real line $0\le \eta \le 2$ ($\eta _0=2)$ with boundary condition $F_1(\eta =2)$. The solution value at $\eta =0$ was computed as part of the overall solution. In figure 17(a) the square root of the $L_2$ error norm (with respect to $F_1(\eta )$) is shown, where it is seen it is seen that the convergence is first order. Also shown in figure 17(b) is the result of an attempted solution of (C7). There is good agreement with $F_2(\eta )$ given by (C11a,b) for $\eta > 0.5$ but oscillatory behaviour with large amplitude for smaller values. This behaviour was found for all $N$ and for different $\eta _0$.