11institutetext: Marco Fasondini 22institutetext: School of Computing and
Mathematical Sciences
University of Leicester
LE1 7RH Leicester
United Kingdom
22email<EMAIL_ADDRESS>33institutetext: J.R. King
44institutetext: School of Mathematical Sciences and Synthetic Biology
Research Centre
University of Nottingham
NG7 2RD Nottingham
United Kingdom
44email<EMAIL_ADDRESS>55institutetext: J.A.C. Weideman
66institutetext: Department of Mathematical Sciences
Stellenbosch University
Stellenbosch 7600
South Africa
66email<EMAIL_ADDRESS>
# Blow up in a periodic semilinear heat equation

Marco Fasondini  John King  J.A.C. Weideman

###### Abstract

Blow up in a one-dimensional semilinear heat equation is studied using a
combination of numerical and analytical tools. The focus is on problems
periodic in the space variable and starting out from a nearly flat, positive
initial condition. Novel results include various asymptotic approximations
that are, in combination, valid over the entire space and time interval right
up to and including the blow-up time. Preliminary results on continuing a
numerical solution beyond the singularity are also presented.

###### Keywords:

Nonlinear blow up, complex singularities, matched asymptotic expansions,
Fourier spectral methods

##  1 Introduction

The blow-up phenomenon in nonlinear evolution equations has been studied
extensively in the literature. Some studies focus on physical applications,
such as singularity formation in fluids braun2 ; braun1 ; hocking ;
Lushnikov21 , runaway in thermal processes Dold91 ; Herrero ; Lacey83 , and
biological applications jabbari2013discrete . Others deal with numerical
aspects such as the computation of blow-up profiles, estimation of blow-up
times, and singularity tracking Berger88 ; Budd96 ; Keller1993 ; Tourigny94 ;
W03 . For a general review, see Galaktionov2002 .

The present paper is a continuation of the papers Keller1993  and W03 . The
equation considered in these papers is the nonlinear heat equation

$u_{t}=u_{xx}+u^{2}$ (1)

(although more general nonlinearities and more than one space dimension were
also considered in Keller1993 ). In this paper we consider the initial
condition

$u(x,0)=\frac{1}{\alpha-\epsilon\cos x},\quad 0<\epsilon\ll\alpha,$ (2)

i.e. periodic, positive, and nearly flat. (A broader class of nearly flat
initial data is considered in Appendix A.) Figure 1 shows a typical blow-up
scenario for this case.

Figure 1: Typical solution profiles of equation (1), in the case of a nearly
flat initial condition such as (2). The solution stays almost flat for a long
while until a point blow-up occurs over a relatively short period of time.

In W03 , it was shown by numerical computation that the approach to blow up is
not necessarily uniform. That is, there may be times when the diffusive term
dominates, leading to a flattening of the solution profile. At other times,
particularly near blow up, the nonlinearity dominates, leading to a steepening
of the profile. By numerically continuing the solution into the complex plane,
it was shown that this behaviour can be associated with the dynamics of the
complex singularities of the solution. When diffusion dominates, they
typically move farther from the real axis while the opposite is true when
nonlinearity dominates. A point blow up occurs when the singularities reach
the real axis. For the nearly flat initial conditions considered in this
paper, however we show that there is no simple correspondence between the
flatness of the solution profile on the real interval $[-\pi,\pi]$ and the
proximity of the nearest singularity, except in the blow-up limit, where the
steepening solution profile is associated with the impingement of
singularities onto the real axis at $x=0$.

In Keller1993 , the focus was not on complex singularities but on asymptotic
estimation of the blow-up time and the solution profile near blow up. The
starting point in that paper was the substitution $u=1/v$, which transforms
(1) into

$v_{t}=v_{xx}-1-2(v_{x})^{2}/v.$ (3)

This transformation has advantages for both analysis and computation.
Advantages for analysis are spelled out in Dold91 . For numerical computation,
it is easier to deal with solutions tending to zero than to infinity. This
avoids the need for specialized rescaling algorithms or moving mesh methods
Berger88 ; Budd96 . The downside is that the simple polynomial nonlinearity in
(1) has been replaced by the more complicated nonlinearity in (3). While this
can be ameliorated by multiplication by $v$, it raises the red flag of
division by zero as $v$ becomes small. However, it was shown in Keller1993
and reaffirmed in section 3 of the present paper that the nonlinear term
remains bounded as $v$ approaches zero. This further suggests the possibility
of integration through the zero of $v$, i.e. through the singularity of $u$,
but Keller1993  reported a failed effort. We continued those investigations
and announce preliminary findings here.

In W03  periodic boundary conditions were considered while Keller1993  looked
at the pure initial value problem on the entire real line. In this paper we
continue with the periodic case, as this gives us access to highly accurate
Fourier spectral methods (which can also be applied to problems on the
infinite line, but with a substantial penalty in accuracy). One contribution
here is a conversion of the analysis of Keller1993  to the periodic situation,
which is not just a triviality but contributes significant new results as we
shall discuss.

For the numerical computations of this paper, a full spectral method based on
a Fourier series

$v(x,t)=\sum_{k=-\infty}^{\infty}c_{k}(t)e^{ikx},\qquad-\pi\leq x<\pi,$ (4)

is used. The derivatives on the right-hand side of (3) are computed by
analytical differentiation of this series, and the nonlinear terms by
convolution and de-convolution. The result is an infinite system of ODEs for
the evolution of the Fourier coefficients $c_{k}$, which we truncate at
$|k|=128$. To integrate this system, we use ode45, MATLAB’s standard ODE
solver. It has adaptive time stepping to maintain accuracy, and the tolerance
parameters for this are set to a stringent $10^{-12}$. To compute the time
$t=t_{c}$ at which the numerical solution blows up (i.e. $v=0$), we use the
fact that the initial conditions considered here lead to blow up at $x=0$ and
hence we check when the sum of the $c_{k}$ equals zero. This can be done by
the ‘event’ option in ode45.

We use this numerical solution as reference solution for the purposes of
checking the various asymptotic estimates. Because of the relatively smooth
nature of $v$, even at the critical time (as will be discussed), we believe it
is sufficiently close to the true solution for all verification purposes.

The three main sections of the paper can be summarized in a nutshell as:
before blow up, at blow up, and after blow up. More specifically, in section 2
we present a perturbation analysis that approximates the solution to (3)
accurately to $O(\epsilon^{2})$ on the entire periodic space domain and the
whole time interval until a time that is $O(\epsilon)$ close to blow up.
Beyond that time up to the blow-up point it has to be modified and this is
done by matched asymptotic expansions, the details of which are contained in
Appendix A. In this appendix we also analyse the dynamics of the singularities
of the solution and in Appendix B the relation between the proximity of the
singularities to the real axis and the steepness of the solution profile on
$[-\pi,\pi]$ is clarified. In sections 2 and 3, the analyses in the appendices
are confirmed by numerical experiments. Section 4 is a preliminary report on
our efforts in integrating through the singularity at the blow-up time and the
subsequent evolution.

##  2 Two-mode perturbation analysis

The analysis of Keller1993  was based on the truncated Taylor expansion

$v\approx a(t)+b(t)x^{2}.$ (5)

By substituting into (3) and dropping powers of $x^{4}$ the problem was
reduced to a dynamical system in the variables $a$ and $b$; see (18) below.
Here we follow an analogous procedure, but consider instead a truncated
Fourier expansion

$v\approx a(t)-b(t)\cos x.$ (6)

Under the assumption of strictly positive solutions, i.e., $0<b(t)<a(t)$, both
of these approximations blow up (in the variable $u=1/v$) in finite time at
$x=0$.

Substitution of (6) into (3) and neglecting $\cos 2x$ contributions gives the
system

$a\frac{da}{dt}+\frac{1}{2}b\frac{db}{dt}=-a-\frac{3}{2}b^{2},\quad
b\frac{da}{dt}+a\frac{db}{dt}=-ab-b,$ (7)

or, in explicit form,

$\frac{da}{dt}=\frac{2ab^{2}+2a^{2}-b^{2}}{b^{2}-2a^{2}},\quad\frac{db}{dt}=\frac{b(2a^{2}-3b^{2})}{b^{2}-2a^{2}}.$
(8)

(Note that the assumption $0<b<a$ precludes the vanishing of the
denominators.) The phase plane of this system is shown in Figure 2.

Figure 2: Phase plane of the system (8) in the domain $0<b<a$. The analysis of
this section approximates solution curves close to $b=0$. Blow up at $x=0$ in
the $u$-equation (1) corresponds to solution trajectories intersecting the
dashed line $b=a$.

Consider solution curves in Figure 2 that originate near $b=0$, say

$a(0)=\alpha,\quad b(0)=\epsilon,\quad 0<\epsilon\ll\alpha.$ (9)

As an explicit solution of the system (7) appears not to exist, we settle for
a perturbation analysis, by expanding

$a=a_{0}(t)+\epsilon a_{1}(t)+O(\epsilon^{2}),\qquad b=\epsilon
b_{1}(t)+O(\epsilon^{2})$ (10)

where $a_{0}(0)=\alpha$, $b_{1}(0)=1$, $a_{1}(0)=0$. Substitution into (7)
gives, to zeroth order

$a_{0}\frac{da_{0}}{dt}=-a_{0}\quad\Rightarrow\quad a_{0}=\alpha-t.$ (11)

At $O(\epsilon)$,

$\frac{da_{1}}{dt}=0,\ \frac{db_{1}}{dt}=-b_{1}\quad\Rightarrow\quad a_{1}=0,\
b_{1}=e^{-t}.$ (12)

This gives, to $O(\epsilon^{2})$, the approximate solution

$\widetilde{v}=\alpha-t-\epsilon\,e^{-t}\cos x.$ (13)

Plugging this expression into (3) gives

$\widetilde{v}\big{(}\widetilde{v}_{t}-\widetilde{v}_{xx}+1+2(\widetilde{v}_{x})^{2}/\widetilde{v}\big{)}=2\,\epsilon^{2}e^{-2t}\sin^{2}x,$
(14)

which confirms that the $v$-equation is satisfied to $O(\epsilon^{2})$
uniformly in $x$, for all $\widetilde{v}$ bounded away from zero. Figure 3
illustrates that the perturbation approximation (13) and the (numerically
computed) solution to the two-mode system (8) are good approximations to the
numerical reference solution of (3) on most of the interval $[0,t_{c})$, where
$t_{c}$ is the critical time at which blow up occurs. As $t\to t_{c}$, the
assumption underlying the perturbation analysis and the two-mode
approximation, $b(t)\ll a(t)$, is no longer valid, however, and the
approximations lose accuracy.

Figure 3: The maximum relative error on $[-\pi,\pi]$ of the perturbation
approximation (13) (blue) and the numerical solution to the two-mode system
(8) (red) as a function of $t$. The errors are calculated with reference to
the numerical solution of (3) mentioned in the final paragraph of section 1.

Various estimates can be obtained from the perturbation solution (13). By
setting $\widetilde{v}=0$ at $x=0$ and excluding $O(\epsilon^{2})$ terms, for
example, one obtains the following estimate for the blow-up time

$t_{c}\approx\widehat{t}_{c}:=\alpha-e^{-\alpha}\epsilon.$ (15)

In Appendix A, using the method of matched asymptotic expansions, a higher-
order estimate of $t_{c}$ is derived (see (43) and (67)):

$t_{c}\approx\widetilde{t}_{c}:=\alpha-e^{-\alpha}\epsilon-(2C_{1}+C_{2}+C_{3})\epsilon^{2},$
(16)

the constants $C_{i}$ being defined in (66). The accuracy of these estimates
is verified in Table 1. (The values of $t_{c}$ listed in the table were
computed by the method described in section 1, and are believed to be correct
to all digits shown.)

| $\alpha=0.25$
---|---
$\epsilon$ | $t_{c}$ | $t_{c}^{\prime}-t_{c}$ | $\widehat{t}_{c}-t_{c}$ | $\widetilde{t}_{c}-t_{c}$
0.1 | 0.161963 |  $-$3.6e-04 | 1.0e-02 | 2.6e-02
0.01 | 0.242093 | 1.8e-03 | 1.2e-04 | 2.8e-04
0.001 | 0.249220 | 2.2e-04 | 1.2e-06 | 2.8e-06
| $\alpha=1$
0.1 | 0.955542 | 2.1e-03 | 7.7e-03 | 4.5e-03
0.01 | 0.996241 | 9.5e-04 | 8.1e-05 | 4.9e-05
0.001 | 0.999631 | 1.1e-04 | 8.1e-07 | 4.9e-07
| $\alpha=4$
0.1 | 3.996685 | 5.0e-04 | 1.5e-03 | 1.4e-05
0.01 | 3.999802 | 5.3e-05 | 1.5e-05 | 1.2e-07
0.001 | 3.999982 | 5.4e-06 | 1.5e-07 | 1.2e-09
Table 1: Blow-up times for various parameter choices in the initial condition
(2): $t_{c}$ is the blow-up time as computed from the reference solution,
$t_{c}^{\prime}$ is the blow-up time as estimated from a numerical solution of
the two-mode system (7) and $\widehat{t}_{c}$ and $\widetilde{t}_{c}$ denote
the estimates (15) and (16), respectively.

The singularity dynamics mentioned in the introduction can be estimated as
follows. The complex singularity in the $u$-equation corresponds to a complex
zero of the $v$-equation. Assuming it is located at $x=iy$ with $y$ real,
setting $\widetilde{v}=0$ gives

$y=\cosh^{-1}\big{(}(\alpha-t)e^{t}/\epsilon\big{)}.$ (17)

This approximation is valid, however, only for small values of $y$, i.e., near
blow up. This follows from the fact that the remainder term on the right-hand
side of (14) grows exponentially with $y$. In Figure 4, the asymptotic
estimate (17) is compared to a numerical estimate of the singularity location
obtained via the method of sulem1983 , which is essentially an estimation of
the width of the strip of analyticity of the solution, by examining the rate
of decay of its Fourier coefficients. To apply the method of sulem1983 , we
use the fact that, to leading order and away from the blow-up time,
singularities of solutions to the $u$-equation are second-order poles111In FKW
we show that these singularities are in fact logarithmic branch points,
however the branch point singularity appears in the fourth-order term in the
local expansion about the singularity. Therefore, for the purpose of
estimating the position of the singularity, we only use its leading-order,
second-order pole behaviour. We shall find in section 3 that in the limit
$t\to t_{c}$, the leading order behaviour of the singularity at $x=0$, which
results from the coalescence of two singularities, is of a more complicated
form than that of a second-order pole. . Figure 4 also shows asymptotic
estimates of the singularity location that are derived in Appendix A.4 by the
method of matched asymptotic expansions. The estimates are valid in the limits
$t\to 0^{+}$ and $t\to t_{c}^{-}$ but in between, for $t=O(1)$, the
singularity location of the original PDE (3) is described by the singularity
location of a more difficult nonlinear initial value problem (68)–(70) which
is not analytically solvable. The asymptotics nonetheless correctly indicate
the initial movement of the singularity away from the real axis (at a speed
that becomes infinite as $t\to 0^{+}$) and the final motion shortly before the
singularity collides with the real axis at $t=t_{c}$.

Figure 4: The position on the positive imaginary axis of the singularity of
the solution to (3) with initial condition $v(x,0)=\alpha-\epsilon\cos x$. (In
the right frame, the approximation (72) is not visible because it is
indistinguishable from (74).)

In Appendix B we investigate the relationship between the position of the
complex singularity and the height of the peak of the solution profile in the
$u$-variable (which is located at $x=0$, see Figure 1) relative to the
solution value at $x=\pm\pi$.

Returning to the analysis of Keller1993 , which was based on the truncated
Taylor approximation (5), we note that the focus in that paper was on the
behaviour at $t=t_{c}$, not the evolution on $[0,t_{c})$ as is the focus here.
It is therefore instructive to adapt that analysis here and compare results to
(13) and (14).

The dynamical system analogous to (7)–(8) is now

$\frac{da}{dt}=2b-1,\quad a\frac{db}{dt}=-8b^{2}.$ (18)

This is a much simpler system, and in fact admits a first integral $2\log
b+1/b+b\log a=\mbox{constant}$, although we shall not make use of this result
and neither was it used in Keller1993 . (No such first integral could be found
for (7)–(8).)

Proceeding with a perturbation analysis based on (9)–(10) give
$a_{0}=\alpha-t$, $a_{1}(t)=2$, and $b_{1}(t)=1$. Therefore, excluding
$O(\epsilon^{2})$ terms, $v$ is approximately

$\widetilde{v}=\alpha-t+\epsilon(2t+x^{2})$ (19)

in analogy with (13). Plugging into (3) gives

$\widetilde{v}\big{(}\widetilde{v}_{t}-\widetilde{v}_{xx}+1+2(\widetilde{v}_{x})^{2}/\widetilde{v}\big{)}=8\,\epsilon^{2}x^{2}.$
(20)

Comparing with (14) shows an advantage of the periodic analysis, namely that
its right-hand side is $O(\epsilon^{2})$ uniformly in $x$, whereas the right-
hand side of (20) is $O(\epsilon^{2})$ only if $x^{2}=O(1)$. That is, with
$\widetilde{v}$ bounded away from zero (13) provides a valid approximation
over the entire (periodic) domain, while (19) is valid only near $x=0$.

On the other hand, the advantage of the analysis of Keller1993  is that it
gives a valid description of the structure of the solution close to blow up,
as discussed in Appendix A.

##  3 Solution in the blow-up limit

Appendix A is devoted to an asymptotic analysis of the solution on three time
scales that are progressively closer to the blow-up time, namely $t=O(1)$,
$T=(t-t_{c})/\epsilon=O(1)$ and
$\tau=-\epsilon\log(-T)=-\epsilon\log((t_{c}-t)/\epsilon)=O(1)$. The analysis
is performed using the method of matched asymptotics and the results on the
first time scale are the same as the those obtained in section 2 via a regular
perturbation analysis (in particular, recall the approximation (13) and its
failure close to the blow up time, as seen in Figure 3).

On the second time scale, $T=O(1)$, the following asymptotic approximation is
derived for $x=O(1)$:

$\begin{split}&v\sim
t_{c}-t+2\epsilon\,e^{-\alpha}\sin^{2}(x/2)+2\epsilon^{2}\log\epsilon\,e^{-2\alpha}\sin^{2}x+\epsilon(t-t_{c})e^{-\alpha}\cos
x\\\
&+2\epsilon^{2}\sin^{2}x\bigg{(}e^{-2\alpha}\log\left(\frac{t_{c}-t}{\epsilon}+2e^{-\alpha}\sin^{2}(x/2)\right)+C_{1}+C_{3}\bigg{)},\end{split}$
(21)

which is obtained by combining (44)–(46) and (65)–(67). It also follows that,
by setting $t=t_{c}$ in (21), we obtain a representation of the solution
profile at the blow-up time that is valid as $t\to t_{c}^{-}$, $\epsilon\to 0$
with $x=O(1)$:

$v\sim
2\epsilon\,e^{-\alpha}\sin^{2}(x/2)+2\epsilon^{2}\sin^{2}x\Big{(}e^{-2\alpha}\log\left(2\epsilon\,e^{-\alpha}\sin^{2}(x/2)\right)+C_{1}+C_{3}\Big{)}.$
(22)

Figure 5 verifies the accuracy of the asymptotic approximations (13) (away
from the blow-up time) and (21)–(22) (close to and at the blow-up time). The
left frame of Figure 6 shows the numerical solution at the blow-up time with
the blow-up profile (22) superimposed on it. Hence, we have accurate
asymptotic expressions for the solution on the entire spatial interval
$x\in[-\pi,\pi]$ and from $t=0$ all the way up to and including the blow-up
time. For $x\to 0$ at $t=t_{c}$, however, we shall need another asymptotic
approximation, namely (24), to be discussed below.

Figure 5: Errors in the asymptotic approximations (13) (blue) and (21) (red)
when compared to the reference solution. The numerical integration method
returns approximations at times $t=t_{k}$, $k=0,1,\ldots,M$, where $t_{0}=0$
and $t_{M}=t_{c}$. Here $M=1155$ and $t_{c}=0.9996$ (more digits are listed in
Table 1). At each $t=t_{k}$, the relative errors of (13) and (21) are
calculated on the interval $x\in[-\pi,\pi]$ and plotted against the time step
index $k$. Because of the adaptive time-stepping, the $t_{k}$ values are not
equidistant, but the spacing is much denser near $t=t_{c}$. As an indication
of this, note that the two curves intersect at $k=372$ and $t_{372}\approx
0.973$, which is already quite close to the critical time $t_{c}$ even though
roughly 800 more time steps are to be taken. From the asymptotic analysis, we
expect the approximation on the first time scale (13) (blue) to start to break
down and (21) (red) to be valid on the second time scale when
$(t_{c}-t)/\epsilon=O(1)$. Here we have that at around $k=800$,
$t_{c}-t_{k}\approx\epsilon=0.001$, which is indeed when the error curves of
the two approximations start to diverge noticeably.

The rate of decay of the Fourier coefficients of $v$ will be of relevance in
the next section when we consider the possibility of continuing the solution
beyond blow-up. From (22) we deduce that the $c_{k}$ decay as $O(k^{-3})$
close to the blow-up time because

$\frac{1}{2\pi}\int_{-\pi}^{\pi}\sin^{2}\\!x\,\log\left(\sin^{2}(x/2)\right)e^{-ikx}dx=\frac{2}{k(k^{2}-4)}$

which implies that

$c_{k}(t_{c})\sim\frac{4\epsilon^{2}e^{-2\alpha}}{k^{3}},\qquad k\to\infty.$
(23)

The right frame of Figure 6 confirms the accuracy of this estimate as
$\epsilon\to 0$.

Figure 6:  Left: The solution to the $v$-equation (3), displayed on a log-
scale, at the blow-up time $t_{c}$ (solid curve) compared to the asymptotic
estimate (22) (dashed curve on top of the solid curve). Here we would like to
draw attention to the fact that, at the origin and at the critical time, the
values of $v$ are computed to levels close to roundoff error ($\sim
10^{-16}$), which translates into huge values in the $u$ solution near blow
up. (See for example the third frame in Figure 9.) This is a consequence of
solving equation (3) instead of (1). Right: The Fourier coefficients of the
reference numerical solution at the blow-up time (dots) compared to the
coefficients of the global asymptotic blow-up profile (dashed lines), given in
(23). The single solid line shows the asymptotic estimate of the decay of the
Fourier coefficients of the local blow-up profile, which is given in (25).

At the blow-up time, for $x$ exponentially small with respect to $\epsilon$,
(22) is no longer valid (see the discussion below (61)). Instead,

$v\sim\frac{\epsilon\,e^{-\alpha}x^{2}}{2-8\epsilon\,e^{-\alpha}\log(x^{2})},$
(24)

which follows from (62) and (67). Figure 7 confirms that (24) has better
accuracy than (22) for small222Ideally, we would show the accuracy of (24) for
even smaller values of $x$ than those in Figure 7. However, this would require
high-precision computations. In standard double precision (with a machine
precision of approximately $10^{-16}$) we cannot compute the errors for
smaller $x$ because round-off errors prevent the accurate evaluation of $v$
from its numerically computed Fourier coefficients. $x$ at the blow-up time.

Figure 7: The relative error of the asymptotic approximations (22) (blue) and
(24) (red) at the blow-up time for small $x$.

Regarding the strength of the singularity at $x=0$, it follows from (24) that
$v^{(n)}\sim 0$, $x\to 0$ for $n=0,1,2$ but the third derivative, however,
blows up as $v^{(3)}=O\left(\frac{1}{x\log^{2}|x|}\right)$, $x\to 0$. A
singularity with two bounded derivatives is consistent with the $O(k^{-3})$
decay of the Fourier coefficients that was derived from (22)333The ‘global’
blow-up profile (22) has a stronger singularity at $x=0$ than the ‘local’
blow-up profile (24) since its second derivative blows up logarithmically at
$x=0$ while for the local profile $v^{(2)}\sim 0$, $x\to 0$. The Fourier
coefficients of the ‘global’ profile nevertheless decay at the correct
$O(k^{-3})$ rate since its second derivative, though unbounded, is
integrable.. In fact, the asymptotic approximation (24) suggests that the
Fourier coefficients of the solution precisely at the blow-up time decay
slightly faster than the $O(k^{-3})$ suggested by (23). To show this, symmetry
and integration by parts can be used to obtain

$c_{k}(t_{c})=\frac{1}{\pi}\int_{0}^{\pi}v\cos kx=\frac{1}{\pi
k^{3}}\int_{0}^{\pi}v^{\prime\prime\prime}\sin kx\,dx.$

From (24) it follows that $v\sim-x^{2}/(8\log(x^{2}))$,
$v^{\prime\prime\prime}\sim 1/(8x\log^{2}(x))$, and by making the change of
variable $x=\mu/k$ with $k\gg\mu$, one obtains,

$c_{k}(t_{c})\sim\frac{1}{8\pi
k^{3}\log^{2}k}\int_{0}^{\infty}\frac{\sin\mu}{\mu}d\mu=\frac{1}{16k^{3}\log^{2}k}.$
(25)

The solid line in Figure 7 shows the estimate (25), which, unlike the estimate
(23), is independent of $\epsilon$ and $\alpha$. The numerical Fourier
coefficients shown in Figure 7 do not decay as fast as (25), which suggests
that one would need to perform high-precision computations (to compute the
solution closer to the blow-up time or with a larger number of Fourier
coefficients) to observe the rate of decay predicted by (25).

For the possible continuation beyond blow-up one needs to confirm that the
nonlinear term in the $v$-equation (3), i.e., $(v_{x})^{2}/v$, remains bounded
at $x=0$ at the blow-up time. In fact, from (24) it follows that the nonlinear
term vanishes, according to $O(1/\log|x|)$, for $x\to 0$.

##  4 Integrating through the singularity

The fact that the nonlinear term in (3) remains bounded as $v$ approaches zero
raises the intriguing possibility of numerically integrating through the blow
up. This was tried in Keller1993 , but the authors found “…the calculation
actually continues the solution slightly beyond the blow-up time of the
original solution $u$. The method becomes unstable a short time after the
blow-up happens, however.” We report here on renewed efforts in this
direction. (Details of how our numerical strategies differ from those of
Keller1993  are postponed to the end of the section.)

Figure 8 shows a series of snapshots of the solution $v$ at different times.
First, observe that there is no sign of instability as the solution passes
smoothly through $v=0$ (third frame). What happens next might be unexpected,
namely, the solution turns complex. In addition, uniqueness is lost: the
complex conjugate of the solution shown is equally likely to appear when
parameters (such as the error tolerances in the time integration scheme or the
number of terms retained in the Fourier series (4)) are adjusted slightly.
Surprising as these results may be, both the non-uniqueness and the fact that
the solution turns complex are consistent with theoretical results of Masuda84
.

Figure 8: Solution to the $v$-equation (3) corresponding to the initial
condition (2) with $\alpha=0.25$, $\epsilon=0.1$. (The corresponding $u$
solution is shown in Figure 9, and a movie of the dynamics can be seen in the
supplementary material that accompanies this paper.) The two noteworthy times
are $t=t_{c}$ (third frame) and $t=2t_{c}$ (seventh frame), where
$t_{c}=0.162$, approximately. At $t=t_{c}$ there is a zero of $v$ at $x=0$,
after which the solution turns complex: the real part is shown in blue and the
imaginary part in red. At $t=2t_{c}$ there are approximate zeros of $v$ near
$x=\pm\pi$. The solution shown is not unique (its conjugate is equally
probable, as discussed in the text.)

Figure 9 shows the $u$-solution that corresponds to the $v$-solution in Figure
8. After the solution turns complex at the critical time (third frame) the
modulus of $u$ shows wave-like behaviour, with two waves travelling in
opposite directions from the origin until they reach the edge of the domain
$x=\pm\pi$. Because of periodicity they meet up with similar waves from
adjacent intervals and a second blow up almost occurs, this time at $x=\pm\pi$
(seventh frame). The modulus grows considerably and we conjecture that by
varying the parameters in the initial condition a proper secondary blow up may
be found. Computing over a longer time interval suggests that $u$
asymptotically approaches the constant (and real) solution $u\sim-1/t$ as
$t\to\infty$.

Figure 9: Same as Figure 8 but here the $u$-solution is shown. Note the blow
up at $t=t_{c}$ and the near blow up at $t=2t_{c}$. The dashed curve is the
modulus $|u|$ (shown only after the first blow up). Note that the scales on
the vertical axes are different in each frame.

We present the results of Figures 8 and 9 knowing full well that we have
little theory to draw on as validation. The theoretical results of Masuda84
suggest the possibility of a complex and non-unique solution post blow-up, but
does not allow for a quantitative comparison. Nevertheless, the following
heuristic observations support the validity of the results shown here.

Firstly, the fact that the numerical solution approaches $-1/t$ as
$t\to\infty$ provides some confidence, as this solves (1) exactly. Secondly,
as noted before, the singularity in the $v$-equation at the critical time is
rather weak, leading to Fourier coefficients that decay at the relatively
rapid rate $c_{k}=O(1/(k^{3}\log^{2}k))$; recall (25). While still much slower
than the typical exponential decay rate for analytic periodic functions, we
conjecture that this decay is nevertheless sufficiently rapid that the
accuracy loss of the spectral method at the critical time is not disastrous.
Figure 10 shows the Fourier coefficients near and at the critical time.

As far as we know, the only published results that deal with numerical
computations of post-blow-up solutions are Cho16 ; Takayasu22 . These authors
based their computations on complexification of the $t$-variable. Following
this idea we integrated along a path that contains a semi-circle in the
complex $t$-plane, centred at the estimated singularity. Using this approach
we obtained the same results as those shown in Figures 8 and 9.

It should be noted that the methods of Cho16 ; Takayasu22  are based on the
$u$-equation, not the $v$-equation. Because the solution to the $u$ equation
grows without bound near the critical time, these methods cannot compute blow-
up solutions in the immediate neighbourhood of the critical time accurately.
For the same reason, estimates of blow-up times are unreliable. Adapting the
methods of Cho16 ; Takayasu22  from the $u$-equation to the $v$-equation might
be a worthwhile future project.

Figure 10: Modulus of the Fourier coefficients $c_{k}$ of the solution $v$
shown in Figure 8. The times are just prior to $t=t_{c}$, at $t=t_{c}$, and
shortly after. (The value of $t_{c}$ is given to more digits in Table 1.) In
the first frame the coefficients decay exponentially, indicative of a function
analytic in a neighbourhood of the real axis. At $t=t_{c}$ the coefficients
decay algebraically as given by (25), because of the influence of the
singularity of $v$ as it reaches the origin. Immediately afterwards the
exponential decay is recovered.

The observation of a complex solution after the critical time has lead us to
conjecture that the failure of the finite difference method of Keller1993  was
because it was almost surely coded as a real system using real arithmetic. One
way to allow for possible complex solutions is to use complex arithmetic
(default in MATLAB) or simply to split the $v$-equation (3) into its real and
imaginary parts. We integrated such a separated system with the Fourier
spectral method mentioned in the first section, which is how the results of
Figures 8–10 were computed. The same idea applies, however, to the finite
difference method of Keller1993 .

Initialising the imaginary part to strictly zero values, the ODE software
signals a singularity at the critical time. Initializing it with a small,
random perturbation on the order of machine roundoff level ($\sim 10^{-16}$),
however, allowed the solution $v$ to pass through the critical $v=0$ without
any sign of instability, as shown in Figure 8.

One priority for future investigations is a better understanding of the
transition from a real to a complex solution and the associated non-
uniqueness. In our case the nonzero imaginary part is triggered by noise at
the level of roundoff error. The radomness dictates whether the continuation
is with one solution or with its complex conjugate.

##  5 Conclusions

We investigated, asymptotically and numerically, point blow-up solutions to a
periodic nonlinear heat equation (1) with nearly flat initial data by
considering the solution $u$ in the reciprocal variable $v=1/u$. We derived
asymptotic approximations for the solution on the entire spatial interval and
from $t=0$ up to and including at the blow-up time, for which we also derived
a second-order approximation. Due to the high accuracy of the Fourier spectral
method (including at the blow-up time at which the $v$-solution has a weak
singularity, unlike the $u$-solution), we were able to check numerically the
validity of the asymptotics. We believe it is unusual for a single numerical
method to confirm asymptotics in multiple regimes since typically numerical
methods are weak, or highly inefficient, in most asymptotic limits unless they
are highly specialised. The key to the success of the numerical method used
here is the fact that it approximates the relatively well-behaved $v$-equation
(3) rather than the $u$-equation (1) whose solution becomes unbounded.

The investigations in this paper point to a number of topics for future
research, not least of which is the validity of the post-blow-up solutions
computed in section 4. In addition, the dynamics of complex singularities of
blow-up solutions to (1) for a larger class of initial conditions, including
initial data leading to non-generic forms of blow up, will be investigated in
FKW , also via a combination of asymptotic and numerical methods. In FKW  we
shall also explore the singularity structure of these blow-up solutions on
their Riemann surfaces in the complex $x$-plane.

###### Acknowledgements.

The first author is grateful to Saleh Tanveer for stimulating and insightful
discussions. This research was started while the three authors were in
residence at the Isaac Newton Institute for Mathematical Sciences as part of
the programme Complex analysis: techniques, applications and computations.
This programme was supported by: EPSRC grant number EP/R014604/1. The work of
the first author was also supported by the Leverhulme Trust Research Project
Grant RPG-2019-144. A grant to the third author from the H.B. Thom foundation
of Stellenbosch University is also gratefully acknowledged.

## References

  * [1] M. Berger and R. V. Kohn.  A rescaling algorithm for the numerical calculation of blowing-up solutions.  Comm. Pure Appl. Math., 41(6):841–863, 1988.
  * [2] S. Braun and A. Kluwick.  Unsteady three-dimensional marginal separation caused by surface-mounted obstacles and/or local suction.  J. Fluid Mech., 514:121–152, 2004.
  * [3] S. Braun and A. Kluwick.  Blow-up and control of marginally separated boundary layers.  Philos. Trans. R. Soc. Lond. Ser. A Math. Phys. Eng. Sci., 363(1830):1057–1067, 2005.
  * [4] C. J. Budd, Jianping Chen, W. Huang, and R. D. Russell.  Moving mesh methods with applications to blow-up problems for PDEs.  In Numerical analysis 1995 (Dundee, 1995), volume 344 of Pitman Res. Notes Math. Ser., pages 1–18. Longman, Harlow, 1996.
  * [5] C.-H. Cho, H. Okamoto, and M. Shoji.  A blow-up problem for a nonlinear heat equation in the complex plane of time.  Jpn. J. Ind. Appl. Math., 33(1):145–166, 2016.
  * [6] J. W. Dold.  On asymptotic forms of reactive-diffusive runaway.  Proc. Roy. Soc. London Ser. A, 433(1889):521–545, 1991.
  * [7] M. Fasondini, J. King, and J.A.C. Weideman.  Complex-plane singularity dynamics for blow up in a nonlinear heat equation: analysis and computation.  In preparation.
  * [8] V. A. Galaktionov and J. J. L. Velázquez.  The problem of blow-up in nonlinear parabolic equations.  Disc. Cont. Dyn. Syst., 8(2):399–433, 2002.
  * [9] M. A. Herrero and J. J. L. Velázquez.  Plane structures in thermal runaway.  Israel J. Math., 81(3):321–341, 1993.
  * [10] L.M. Hocking, K. Stewartson, J.T. Stuart, and S.N. Brown.  A nonlinear instability burst in plane parallel flow.  J. Fluid Mech., 51(4):705–735, 1972.
  * [11] S. Jabbari and J.R. King.  Discrete and continuum multiscale behaviour in bacterial communication.  In Multiscale Computer Modeling in Biomechanics and Biomedical Engineering, pages 299–320. Springer, 2013.
  * [12] J. B. Keller and J. S. Lowengrub.  Asymptotic and numerical results for blowing-up solutions to semilinear heat equations.  In Singularities in fluids, plasmas and optics (Heraklion, 1992), volume 404 of NATO Adv. Sci. Inst. Ser. C Math. Phys. Sci., pages 111–129. Kluwer Acad. Publ., Dordrecht, 1993.
  * [13] A. A. Lacey.  Mathematical analysis of thermal runaway for spatially inhomogeneous reactions.  SIAM J. Appl. Math., 43(6):1350–1366, 1983.
  * [14] P. M. Lushnikov, D. A. Silantyev, and M. Siegel.  Collapse versus blow-up and global existence in the generalized Constantin-Lax-Majda equation.  J. Nonlinear Sci., 31(5):Paper No. 82, 56, 2021.
  * [15] K. Masuda.  Analytic solutions of some nonlinear diffusion equations.  Math. Z., 187(1):61–73, 1984.
  * [16] C. Sulem, P.L. Sulem, and H. Frisch.  Tracing complex singularities with spectral methods.  J. Comput. Phys., 50(1):138–161, 1983.
  * [17] A. Takayasu, J.P. Lessard, J. Jaquette, and H. Okamoto.  Rigorous numerics for nonlinear heat equations in the complex plane of time.  Numer. Math., 151(3):693–750, 2022.
  * [18] Y. Tourigny and M. Grinfeld.  Deciphering singularities by discrete methods.  Math. Comp., 62(205):155–169, 1994.
  * [19] J.A.C. Weideman.  Computing the dynamics of complex singularities of nonlinear PDEs.  SIAM J. Appl. Dyn. Syst., 2:171–186, 2003.

##  Appendix A Analysis via the method of matched asymptotic expansions

###  A.1 Truncated Fourier expansion

In this section we revisit and expand the analysis based on the two-mode
Fourier truncation outlined in section 2. Recall that by substituting (6) into
the $v$-equation (3) and neglecting the $\cos 2x$ term the system (7) was
obtained. From this, one finds by a self consistency argument that near blow-
up, generically (i.e. even for $\epsilon=O(1)$)

$a\sim a_{c}+(1+2a_{c})(t_{c}-t),\qquad b\sim a_{c}-a_{c}(t_{c}-t),\qquad t\to
t_{c}^{-},$ (26)

where $a_{c}$ and the blow up time $t_{c}$ are positive constants, so that the
blow-up behaviour associated with (6) takes the form

$v\sim(1+3a_{c})(t_{c}-t)+\frac{1}{2}a_{c}x^{2},\qquad t\to t_{c}^{-},\qquad
x=O\left((t_{c}-t)^{1/2}\right).$ (27)

We emphasise that (26)–(27) describe the blow-up behaviour of (7) rather than
that of (3): (27) has features in common with, but is not a valid
representation of, the blow-up behaviour of the full PDE (3), a key point to
which we shall return.

In keeping with our goal of characterising the behaviour of (3) for near-flat
initial conditions, we now develop a fully analytic asymptotic description of
the behaviour of the two-mode system (7) for initial conditions (9) with
$\alpha=O(1)$. There are two timescales; the first coincides with (10), so
that at leading order (11) and (12) follow.

On the second time scale, and near blow up, we set

$t=\alpha+\epsilon T,\quad b=\epsilon B,\quad a=\epsilon
A_{0}(T)+O(\epsilon^{2}),\quad B=B_{0}(T)+O(\epsilon),$ (28)

with $T=O(1)$, so that

$\begin{split}A_{0}\frac{dA_{0}}{dT}+\frac{1}{2}B_{0}\frac{dB_{0}}{dT}&=-A_{0},\\\
B_{0}\frac{dA_{0}}{dT}+A_{0}\frac{dB_{0}}{dT}&=-B_{0}.\end{split}$ (29)

Since $A_{0}\pm\frac{1}{\sqrt{2}}B_{0}$ each satisfy

$\Phi\frac{d\Phi}{dT}=-\Phi$

we find on matching with (11)–(12) that

$A_{0}=-T,\qquad B_{0}=e^{-\alpha}$ (30)

(a result that relies on the observation that $a$ has no $O(\epsilon)$
contribution for $t=O(1)$, see (12)). Thus

$t_{c}\sim\alpha-\epsilon e^{-\alpha},\qquad a_{c}\sim\epsilon
e^{-\alpha},\qquad\epsilon\to 0^{+}$

in (26)–(27). Moreover, while the rescaling (28) leads to a modified balance,
i.e. (29) in place of (11) and (12), the result (30) implies that the solution
passes unscathed through $T=O(1)$, with (11)–(12) being recovered for
$t>\alpha$. Thus with the two-mode approximation (6) and (8) continuation
through blow up seems to be straightforward, in contrast to that of the
$v$-equation (3), as we shall subsequently demonstrate.

###  A.2 Truncated Taylor expansion

Much of what follows in this subsection revisits results from [12], which we
derive using a different approach (namely, matched asymptotic expansions).

As in section 2, we also consider the truncated Taylor approximation (5),
which gives rise to the system (18) whose blow-up behaviour takes the form
(again, by a self-consistency argument)

$a\sim t_{c}-t,\qquad b\sim\frac{1}{8(-\log(t_{c}-t)+b_{c})},\qquad t\to
t_{c}^{-}$ (31)

for constants $t_{c}$ and $b_{c}$. In contrast to the results of the previous
subsection, (18) _does_ capture the blow-up behaviour of the full PDE (3), a
point to which we shall also come back.

We now return to initial conditions of the form (9), with the scalings (10)
applying for $t=O(1)$, so expanding in the form

$a\sim a_{0}+\epsilon a_{1},\qquad b=\epsilon B,\qquad B\sim B_{0}+\epsilon
B_{1},$

(18) implies

$a_{0}=\alpha-t,\quad B_{0}=1,\quad a_{1}=2t,\quad
B_{1}=8\log\left((\alpha-t)/\alpha\right).$ (32)

Under the scalings (28) the leading-order balances do not change and the
results of relevance below read

$A_{0}=2\alpha-T,\quad B_{0}=1,\quad
B_{1}=8\log\left((2\alpha-T)/\alpha\right),$ (33)

the matching into $a_{1}$ in (32) leading to the $2\alpha$ contributions. The
final scale in this case is then more subtle than those above: we set

$a=(T_{c}-T)\hat{a}(\tau),\quad
B=\hat{b}(\tau),\quad\tau=-\epsilon\log(T_{c}-T),$ (34)

where

$T_{c}\sim 2\alpha,\qquad t_{c}\sim\alpha+2\epsilon\alpha,\qquad\epsilon\to
0^{+}.$ (35)

The introduction of $\tau$ in (34) is associated with the expansion for $B$ in
(33) disordering and (18) becomes

$\begin{split}&\epsilon\frac{d\hat{a}}{d\tau}-\hat{a}=2\epsilon\hat{b}-1,\\\
&\hat{a}\frac{d\hat{b}}{d\tau}=-8\hat{b}^{2},\end{split}$ (36)

so that, on matching with (33),

$\hat{a}_{0}=1,\qquad\hat{b}_{0}=\frac{1}{1+8\tau},\qquad\hat{a}_{1}=-\frac{2}{1+8\tau},$
(37)

so that $b_{c}\sim 1/(8\epsilon)$ in (31).

For reasons that will become apparent below, it is instructive to record the
implications of (37) under the scaling

$x=\epsilon^{1/2}X$ (38)

whereby

$v\sim\epsilon\left(T_{c}-T+\frac{\epsilon}{1+8\tau}\left(X^{2}-2(T_{c}-T)\right)\right);$

this will reappear in the analysis below in describing the local blow-up
behaviour, but also confirms that (5) cannot describe the spatial profile at
blow up (i.e. at $\tau=+\infty$).

###  A.3 More general near-flat initial data

We now turn to the derivation of the blow-up behaviour, subjecting (3), the
full PDE, to the initial data

$v(x,0)=\alpha+\epsilon V(x)$

with, again, $0<\epsilon\ll 1$. It is striking that this limit allows a near
complete analytical description of the transition to blow up through the fully
nonlinear regime. Three time scales are required. On the first time scale:

$t=O(1),\qquad v\sim v_{0}(t)+\epsilon v_{1}(x,t)+\epsilon^{2}v_{2}(x,t).$ (I)

All three terms in this expansion are required for what follows. An immediate
result is that

$v_{0}=\alpha-t$

leading to

$\frac{\partial v_{1}}{\partial t}=\frac{\partial^{2}v_{1}}{\partial
x^{2}},\qquad v_{1}(x,0)=V(x).$ (39)

We shall consider general $V(x)$, but in the case of real-valued
$2\pi$-periodic initial conditions (or zero Neumann boundary conditions on a
finite domain), we can Fourier decompose in the usual way, so that

$v_{1}(x,t)=\sum_{k=-\infty}^{\infty}\ a_{k}e^{ikx-k^{2}t}.$

We note this special case for two reasons – firstly for its relevance to
section A.1 and secondly for the obvious observation that the high-frequency
modes are rapidly decaying, providing additional motivation for the analysis
of sections 2 and A.1 but being also in some respects counter-intuitive, given
that point blow up subsequently ensues.

We define

$\Phi_{1}(x)=v_{1}(x,\alpha)$

so that

$v_{1}(x,t)\sim\Phi_{1}(x)-(\alpha-t)\Phi_{1}^{\prime\prime}(x),\qquad
t\to\alpha^{-}.$ (40)

At next order we have

$\frac{\partial v_{2}}{\partial t}-\frac{\partial^{2}v_{2}}{\partial
x^{2}}=-\frac{2}{\alpha-t}\left(\frac{\partial v_{1}}{\partial
x}\right)^{2},\qquad v_{2}(x,0)=0,$ (41)

so that

$\begin{split}v_{2}(x,t)\sim&2\left(\Phi_{1}^{\prime}(x)\right)^{2}\log(\alpha-t)+\Phi_{2}(x)-2\left[\left(\Phi_{1}^{\prime}\right)^{2}\right]^{\prime\prime}(\alpha-t)\log(\alpha-t)\\\
&+\left(4\left(\Phi_{1}^{\prime\prime}\right)^{2}-\Phi_{2}^{\prime\prime}\right)(\alpha-t),\qquad
t\to\alpha^{-},\end{split}$ (42)

where (42) serves to define $\Phi_{2}(x)$.

On the second time scale:

$v=\epsilon w,\qquad w\sim
w_{0}(x,T)+\epsilon\log(1/\epsilon)w_{1}(x,T)+\epsilon w_{2}(x,T).$ (II)

Here $t=t_{c}(\epsilon)+\epsilon T$, $T=O(1)$ with $t_{c}(0)=\alpha$ and with
correction terms to $t_{c}$ identified below444The notation here differs
somewhat from that of the previous subsections.. We shall also need to
consider the rescaling (38). We have

$w\frac{\partial w}{\partial T}=\epsilon\left(w\frac{\partial^{2}w}{\partial
x^{2}}-2\left(\frac{\partial w}{\partial x}\right)^{2}\right)-w$

so, matching with (40), (42), and defining $\beta_{1}$ and $\beta_{2}$ via

$t_{c}(\epsilon)\sim\alpha+\epsilon\beta_{1}+\epsilon^{2}\beta_{2},\qquad\epsilon\to
0$ (43)

it follows for $x=O(1)$ that

$w_{0}=-T+\Phi_{1}(x)-\beta_{1},$ (44)
$w_{1}=-2\left(\Phi_{1}^{\prime}(x)\right)^{2}$ (45)

and hence

$\frac{\partial w_{2}}{\partial
T}=\Phi_{1}^{\prime\prime}(x)-\frac{2}{-T+\Phi_{1}(x)-\beta_{1}}\left(\Phi_{1}^{\prime}(x)\right)^{2},$

implying

$w_{2}=-(-T)\Phi_{1}^{\prime\prime}(x)+2\left(\Phi_{1}^{\prime}(x)\right)^{2}\log\left(-T+\Phi_{1}(x)-\beta_{1}\right)+\Phi_{2}(x)-\beta_{2}.$
(46)

We take blow up to occur at $x=0$, requiring that
$\Phi_{1}^{\prime}(0)=\Phi_{2}^{\prime}(0)=0$ and take

$\Phi_{1}(x)\sim\beta_{1}+\gamma_{1}x^{2},\qquad\Phi_{2}(x)\sim\beta_{2}+\gamma_{2}x^{2},\qquad
x\to 0$ (47)

for constants $\beta_{1,2}$, $\gamma_{1,2}$ with $\gamma_{1}>0$555The case
$\gamma_{1}=0$ corresponds to non-generic forms of blow up – we shall not
pursue such matters here. and where the requirement that $w=0$ at $x=0$, $T=0$
implies that $\beta_{1,2}$ in (43) are specified by (47), with $\Phi_{1,2}$
being determined by the linear problems (39) and (41).

For $X=O(1)$, we first generate the required matching conditions from (40),
(42) and (47), whereby

$\begin{split}w\sim&-T+\left(\epsilon\gamma_{1}-8\epsilon^{2}\log(1/\epsilon)\gamma_{1}^{2}+8\epsilon^{2}\gamma_{1}^{2}\log(-T)+\epsilon^{2}\gamma_{2}\right)\left(X^{2}-2(-T)\right)\\\
&+16\epsilon^{2}\gamma_{1}^{2}(-T),\end{split}$ (48)

wherein we have retained only the required terms in $X^{0}$ and $X^{2}$ – in
general $v_{1}$ will also lead to contributions of the form

$\epsilon^{3/2}\left(X^{3}-6(-T)X\right),\qquad\epsilon^{2}\left(X^{4}-12(-T)X^{2}+12(-T)^{2}\right)$

in $w$ but these can be neglected for our purposes, being sub-dominant as
$T\to 0^{-}$ with $X=O\left((-T)^{1/2}\right)$.

We have

$w\frac{\partial w}{\partial T}=w\frac{\partial^{2}w}{\partial
X^{2}}-2\left(\frac{\partial w}{\partial X}\right)^{2}-w$ (49)

and from this we find that the relevant terms for $T=O(1)$ simply reproduce
the matching condition (48) obtained from the expansion for $t=O(1)$.
Importantly, the expansion (48) disorders for $\epsilon\log(-T)=O(1)$, as in
section A.2, leading us on to our final time scale, as follows:

$\tau=O(1),\quad w=(-T)g,\quad g\sim 1+\epsilon
g_{1}(\xi,\tau)+\epsilon^{2}\log(1/\epsilon)g_{2}(\xi,\tau)+\epsilon^{2}g_{3}(\xi,\tau)$
(III)

where

$\tau=-\epsilon\log(-T),\qquad\xi=\frac{X}{(-T)^{1/2}}.$ (50)

We emphasise that much of what follows reconstructs known blow-up behaviour
(see the review of [8], for example.) The novelty here lies in the focus on
the near-flat initial data, which allows more detailed characterization of the
full spatial behaviour.

The near self-similar solution ansatz (50) transforms (49) to

$\epsilon g\frac{\partial g}{\partial\tau}-g^{2}+\frac{1}{2}\xi
g\frac{\partial
g}{\partial\xi}=g\frac{\partial^{2}g}{\partial\xi^{2}}-2\left(\frac{\partial
g}{\partial\xi}\right)^{2}-g$

so that

$-g_{1}+\frac{1}{2}\xi\frac{\partial
g_{1}}{\partial\xi}=\frac{\partial^{2}g_{1}}{\partial\xi^{2}},$ (51)

implying that

$g_{1}=\sigma_{1}(\tau)\left(\xi^{2}-2\right);$ (52)

similarly

$g_{2}=\sigma_{2}(\tau)\left(\xi^{2}-2\right)$ (53)

but $g_{3}$ satisfies (following cancellation of a number of terms using (51))

$\frac{d\sigma_{1}}{d\tau}\left(\xi^{2}-2\right)-g_{3}+\frac{1}{2}\xi\frac{dg_{3}}{d\xi}=\frac{d^{2}g_{3}}{d\xi^{2}}-8\sigma_{1}^{2}\xi^{2}.$
(54)

Since we need to preclude exponential growth of $g_{3}$ (i.e. to exclude
contributions with $\log g_{3}\sim\xi^{2}/4$ as $\xi\to\pm\infty$), (54) both
requires that

$g_{3}=-2\frac{d\sigma_{1}}{d\tau}+\sigma_{3}(\tau)\left(\xi^{2}-2\right)$
(55)

and generates the solvability condition

$\frac{d\sigma_{1}}{d\tau}=-8\sigma_{1}^{2};$ (56)

that the $O(\epsilon)$ term in $g$ is only fully determined via (54)
necessitates that the expansion be taken up to $O(\epsilon^{2})$. Matching
with (48) then requires that

$\sigma_{1}=\frac{\gamma_{1}}{1+8\gamma_{1}\tau},\qquad-2\frac{d\sigma_{1}}{d\tau}=\frac{16\gamma_{1}^{2}}{(1+8\gamma_{1}\tau)^{2}},$
(57)

simultaneously confirming matching with the second, fourth and final terms in
(48). The calculation of $\sigma_{2}$ and $\sigma_{3}$ requires solvability
conditions at yet higher orders (which we will not pursue), with (48)
requiring that

$\sigma_{2}(0)=-8\gamma_{1}^{2},\qquad\sigma_{3}(0)=\gamma_{2}.$

By analogy with the subdivision into $x$ and $X$ in (II), we need also to
consider $\eta=O(1)$, where

$\eta=\frac{x}{(-T)^{1/2}},\qquad\eta=\epsilon^{1/2}\xi,$

though a third scale (namely $x=O(1)$) will also require consideration here.
Since

$\epsilon g\frac{\partial g}{\partial\tau}-g^{2}+\frac{1}{2}\eta
g\frac{\partial
g}{\partial\eta}=\epsilon\left(g\frac{\partial^{2}g}{\partial\eta^{2}}-2\left(\frac{\partial
g}{\partial\eta}\right)^{2}\right)-g,$

setting

$g\sim G_{0}(\eta,\tau)+\epsilon\log(1/\epsilon)G_{1}(\eta,\tau)+\epsilon
G_{2}(\eta,\tau)$

gives on matching into $\xi=O(1)$

$G_{0}=1+\sigma_{1}(\tau)\eta^{2},\qquad G_{1}=\sigma_{2}(\tau)\eta^{2}$ (58)

and

$\frac{1}{2}\eta\frac{dG_{2}}{d\eta}-G_{2}=2\sigma_{1}-\frac{d\sigma_{1}}{d\tau}\eta^{2}-\frac{8\sigma_{1}^{2}\eta^{2}}{1+\sigma_{1}\eta^{2}},$

so that, again matching into $\xi=O(1)$,

$G_{2}=-2\sigma_{1}(\tau)+8\sigma_{1}^{2}(\tau)\eta^{2}\log\left(1+\sigma_{1}(\tau)\eta^{2}\right)+\sigma_{3}(\tau)\eta^{2};$
(59)

were $\sigma_{1}$ not a solution to (56), a term in $\log\eta$ would also be
present here, further clarifying the status of (56) as a solvability
condition.

Finally, we can simply set $T=0$ in (44)–(46), noting that $\tau=O(1)$
corresponds to exponentially small $T$, to obtain the profile at blow up for
$x=O(1)$. Thus as $t\to t_{c}^{-}$, $\epsilon\to 0$ with $x=O(1)$ we have

$\begin{split}w\sim&\>\Phi_{1}(x)-\beta_{1}-2\epsilon\log(1/\epsilon)\left(\Phi_{1}^{\prime}(x)\right)^{2}\\\
&+\epsilon\left(2\left(\Phi_{1}^{\prime}(x)\right)^{2}\log(\Phi_{1}(x)-\beta_{1})+\Phi_{2}(x)-\beta_{2}\right).\end{split}$
(60)

The expression (60) has small-$x$ behaviour

$w\sim\gamma_{1}x^{2}-8\epsilon\log(1/\epsilon)\gamma_{1}^{2}x^{2}+\epsilon\left(8\gamma_{1}^{2}x^{2}\log(\gamma_{1}x^{2})+\gamma_{2}x^{2}\right),$
(61)

The expressions (60)–(61) do not apply for exponentially small $x$, however.
Instead, we need to extract from (58)–(59) the terms relevant for large
$\eta$, namely

$g\sim\sigma_{1}\eta^{2}+16\epsilon\sigma_{1}^{2}\eta^{2}\log\eta,$

these being of the same order for $\epsilon\log\eta=O(1)$.

Reconstructing $w$ from these using (57), we have

$w\sim\frac{\gamma_{1}}{1+8\gamma_{1}\tau}x^{2}+\frac{16\epsilon\gamma_{1}^{2}x^{2}}{(1+8\gamma_{1}\tau)^{2}}\log\eta,$

which, using

$\tau=-2\epsilon\log x+2\epsilon\log\eta,$

implies that

$w\sim\frac{\gamma_{1}x^{2}}{1+16\epsilon\gamma_{1}\log(1/x)}$ (62)

describes the profile at blow up for $x$ exponentially small with respect to
$\epsilon$ (having taken the various limits in the appropriate order); (62)
matches with the relevant terms in (61) for larger $x$.

We note that non-analytic (i.e. logarithmic) terms here occur as a matter of
course, in $(t_{c}-t)$ in (57) and in $x$ in (61): these arise constructively
in the current analysis rather than being introduced a priori as part of a
solution ansatz. That the current limit provides a detailed asymptotic
characterisation of the profile (60) at blow up for almost all $x$, not just
close to blow-up point, is also noteworthy.

####  A.3.1 Comparison with truncated expansions

The above systematic asymptotic analysis clarifies the extent of applicability
of the ad-hoc approximations treated in sections A.1 and A.2.

Starting with a comparison with the results of section A.2, we have already
noted that the truncated Taylor expansion correctly captures the blow-up
behaviour, as can be substantiated by the following observations. For
$V(x)=x^{2}$ it follows that

$\begin{split}&v_{1}=x^{2}+2t,\qquad
v_{2}=8x^{2}\log\left((\alpha-t)/\alpha\right)-16(\alpha-t)\log\left((\alpha-t)/\alpha\right)-16t,\\\
&\Phi_{1}=x^{2}+2\alpha,\qquad\Phi_{2}=-8x^{2}\log\alpha-16\alpha,\end{split}$
(63)

which imply (see (47) and (43))

$\beta_{1}=2\alpha,\qquad\gamma_{1}=1,\qquad\beta_{2}=-16\alpha,\qquad\gamma_{2}=-8\log\alpha,\qquad
t_{c}\sim\alpha+2\alpha\epsilon-16\alpha\epsilon^{2}.$ (64)

More importantly, (48), (52), (53), (55) and (58) are all quadratic in $x$ for
general $V(x)$, so the approximation in section A.2 represents an attractor in
that sense. As already noted, it cannot capture the spatial blow-up profile,
the first manifestation of this being represented by the logarithmic term in
(59); see also (61) and (62). It is especially noteworthy that (55) is a
quadratic in $\xi$; if the corresponding analysis is undertaken on (1) rather
than on (3), a $\xi^{4}$ term arises at that order. This represents a hidden
benefit of the $v$ formulation.

Turning now to section A.1, corresponding to $V(x)=-\cos x$, we can both
exemplify the form of the profile (60) at blow up and indicate where the
analysis of section A.1 breaks down in describing the behaviour of the full
PDE. In this case

$\begin{split}&v_{1}=-e^{-t}\cos x,\quad
v_{2}=-\int_{0}^{t}\frac{e^{-2s}}{\alpha-s}ds-e^{-4t}\int_{0}^{t}\frac{e^{2s}}{\alpha-s}ds\cos(2x),\
\\\ &\Phi_{1}=-e^{-\alpha}\cos
x,\quad\Phi_{2}=-C_{1}-C_{2}-\left(C_{1}+C_{3}\right)\cos 2x,\end{split}$ (65)

where

$C_{1}=e^{-2\alpha}\log\alpha,\qquad
C_{2}=\int_{0}^{\alpha}\frac{e^{-2t}-e^{-2\alpha}}{\alpha-t}dt,\qquad
C_{3}=e^{-4\alpha}\int_{0}^{\alpha}\frac{e^{2t}-e^{2\alpha}}{\alpha-t}dt,$
(66)

which imply (see (47))

$\beta_{1}=-e^{-\alpha},\qquad\gamma_{1}=\frac{e^{-\alpha}}{2},\qquad\beta_{2}=-2C_{1}-C_{2}-C_{3},\qquad\gamma_{2}=2(C_{1}+C_{3}).$
(67)

That $v_{2}$ and $\Phi_{2}$ contain $\cos 2x$ contributions is already
indicative of the failure of the truncated Fourier expansion close to the blow
up; nevertheless, $v_{1}$ and $w_{0}$ (in (44)) provide the dominant spatial
dependencies on the relevant scales and retain the form of section A.1 (and
the $O(\epsilon)$ term in (48) can be viewed as being associated with the
latter’s Taylor expansion). The truncation is, however, entirely unable to
reproduce the behavior for $\tau=O(1)$, manifesting nothing reflecting the
solvability argument leading to (56).

Notwithstanding their deficiencies, the simple approximations of sections A.1
and A.2 are instructive both for familiar reasons (namely for the transparency
and simplicity of their analysis) and because they can immediately be
analytically continued to assess their relevance to the evolution of
singularities in the complex plane; see Figure 4 (in which the dashed line is
derived from the approximation of section A.1).

###  A.4 Complex-singularity dynamics

The motion of the singularity necessarily breaks down into the three time
scales of Appendix A.3 and here we adopt the specific initial condition
$v=\alpha-\epsilon\cos x$. We set $x=iy$ and describe the location of the
nearest singularity (corresponding to $v=0$) on the positive imaginary axis.

On the first time scale ($t=O(1)$), we set $y=\log(1/\epsilon)+Y$ in (3) and
have the full balance

$v\frac{\partial v}{\partial t}=-\left(v\frac{\partial^{2}v}{\partial
Y^{2}}-2\left(\frac{\partial v}{\partial Y}\right)^{2}\right)-v,$ (68)

now to be solved as an initial value problem subject to

$\displaystyle\text{at }t=0,\qquad\qquad v=\alpha-\frac{1}{2}e^{Y},$ (69)
$\displaystyle\text{as }Y\to-\infty,\qquad v\sim\alpha-t-\frac{1}{2}e^{Y-t}.$
(70)

This not analytically solvable (and in effect a more difficult problem than
the original PDE with nearly flat initial data), but does establish that to
leading order the singularity location is fixed: $y\sim\log(1/\epsilon)$. This
problem can, however, be solved analytically in the limits $t\to 0^{+}$ and
$t\to\alpha^{-}$; using the techniques of [7], we find

$\displaystyle y$ $\displaystyle\sim$
$\displaystyle\log(2\alpha/\epsilon)+\left(2t\log(1/t)\right)^{1/2},\qquad
t\to 0^{+},$ (71) $\displaystyle y$ $\displaystyle\sim$
$\displaystyle\log(2/\epsilon)+\alpha+\log(\alpha-t),\quad\qquad\hskip
2.84544ptt\to\alpha^{-},$ (72)

so the singularity moves away from the real axis at early times (with
unbounded speed as $t\to 0^{+}$) before reversing.

The results on the second and final time scales follow immediately from the
above real line results. On the second time scale ($t=t_{c}+\epsilon T$,
$T=O(1)$), since (see (44) and (65)–(67))

$w_{0}=-T-e^{-\alpha}\left(\cosh y-1\right)$

the singularity location satisfies

$\cosh y\sim 1+e^{\alpha}(-T)$ (73)

so that

$y\sim\log\left(1+e^{\alpha}(-T)+\left(2e^{\alpha}(-T)+e^{2\alpha}(-T)^{2}\right)^{1/2}\right).$
(74)

If follows from (74) that

$\displaystyle y\sim\log(-2T)+\alpha,\qquad T\to-\infty,$ (75) $\displaystyle
y\sim\left(2e^{\alpha}(-T)\right)^{1/2},\qquad T\to 0^{-},$ (76)

where (75) exhibits the necessary matching with (72).

On the final time scale ($\tau=-\epsilon\log(-T)$, $\tau=O(1)$), since (see
(58))

$G_{0}=1-\sigma_{1}y^{2}/(-T),\qquad\sigma_{1}=\frac{\gamma_{1}}{1+8\gamma_{1}\tau},$

with $\gamma_{1}=1/(2e^{\alpha})$ the singularity location satisfies

$y\sim\left(2e^{\alpha}(-T)\right)^{1/2}\left(1-4\epsilon
e^{-\alpha}\log(-T)\right)^{1/2},$ (77)

which matches with (76) for $\epsilon|\log(-T)|\ll 1$, while

$y\sim\left(8(t_{c}-t)\log(1/(t_{c}-t))\right)^{1/2},\qquad t\to t_{c}^{-},$

accordingly describes impingement onto the real axis.

##  Appendix B Flatness of the solution on the real line and proximity of the
nearest singularity

Note from Figure 1 that the solutions to the $u$-equation (1) that we consider
attain their maximum at $x=0$ and their minimum at $x=\pm\pi$. Therefore, as
an indication of the ‘flatness’ of the solution profile, we consider the
quantity $u(0,t)-u(\pi,t):=f(t)$, which is the relative height of the peak of
the solution on $[-\pi,\pi]$. We shall consider the relation between the
flatness of the solution profile and the distance of the singularities to the
real axis.

Let the $a_{k}(t)$ denote the Fourier coefficients of the solution in the $u$
variable. Since the solution is real (prior to blow up) and even,
$a_{k}=a_{-k}$ and thus the flatness of the solution is given by

$\displaystyle
f(t)=\sum_{k=-\infty}^{\infty}a_{k}(t)-\sum_{k=-\infty}^{\infty}(-1)^{k}a_{k}(t)=4\sum_{k=0}^{\infty}a_{2k+1}(t).$
(78)

If $f^{\prime}<0$, the solution becomes more flat while if $f^{\prime}>0$ the
solution becomes steeper on $x\in[-\pi,\pi]$. Using the residue theorem, it
follows that

$\displaystyle a_{k}(0)$ $\displaystyle=$
$\displaystyle\frac{1}{2\pi}\int_{-\pi}^{\pi}\frac{e^{-ikx}}{\alpha-\epsilon\cos
x}dx$ (79) $\displaystyle=$
$\displaystyle\frac{\left(\frac{\alpha}{\epsilon}+\sqrt{\left(\frac{\alpha}{\epsilon}\right)^{2}-1}\right)^{-|k|}}{\epsilon\sqrt{\left(\frac{\alpha}{\epsilon}\right)^{2}-1}}=\frac{1}{\alpha}\left(\frac{\epsilon}{2\alpha}\right)^{|k|}+O(\epsilon^{|k|+2}),\qquad\epsilon\to
0.$

The coefficients satisfy

$a_{k}^{\prime}=-k^{2}a_{k}+b_{k},$ (80)

where $b_{k}$ is the $k$-th Fourier coefficient of $u^{2}$.

It follows from (79)–(80) that, to leading order, as $\epsilon\to 0$, and away
from the blow-up time,

$a_{0}(t)\sim\frac{1}{\alpha-t},\qquad a_{1}(t)\sim\frac{\epsilon
e^{-t}}{2(\alpha-t)^{2}},$ (81)

therefore

$f^{\prime}(t)\sim
4a_{1}^{\prime}\sim\frac{2\epsilon\,e^{-t}(t-(\alpha-2))}{(\alpha-t)^{3}}.$

Hence, on the real line, the solution switches from flattening to steepening
at $t\sim\alpha-2$ if $\alpha>2$ and the minimal flatness is $f(\alpha-2)\sim
4a_{1}(\alpha-2)\sim\epsilon\,e^{2-\alpha}/2$. If $\alpha<2$, the solution
does not flatten at all but steepens from $t=0$ until the blow-up time. Figure
11 confirms the validity of the flatness approximation $f(t)\sim 4a_{1}(t)$
away from the blow-up time.

Figure 11:  Left: The numerically computed flatness $f(t)=u(0,t)-u(\pi,t)$
(solid blue curve) compared to the approximation $f(t)\sim 4a_{1}$, with
$a_{1}$ given in (81) (red dashed curve). Right: the relative error of the
approximation $f(t)\sim 4a_{1}$, which, as expected, ceases to be valid as
blow up is approached.

In the complex plane, we deduce from the asymptotic approximation (71) that,
regardless of the parameter values (provided $\alpha\ll\epsilon$), the
singularities initially move away from the real axis and turn around at
$t\sim\exp(-1)\approx 0.37$, which is consistent with the numerical results in
Figure 4. This illustrates that there is not a simple correspondence between
the distance of the singularity from the real axis and the flatness of the
solution.

Indeed, it follows from (78) and the rapid decay of the Fourier coefficients
(away from the blow-up time) that the flatness of the solution is determined
by the behaviour of the low-order Fourier coefficient $a_{1}$. On the other
hand, the distance of the singularity from the real axis is determined by the
behaviour of the high-order Fourier modes. Hence, the discrepancy between the
flatness of the solution and the proximity of the singularity is due to the
qualitatively different evolution of the lower and higher order Fourier modes.

For the high-order modes, the diffusion term $-k^{2}a_{k}$ dominates the
nonlinear term $b_{k}$ for small $t$ (again, due to the rapid decay of the
Fourier coefficients). Therefore $a_{k}^{\prime}<0$, the high-order Fourier
coefficients decay for small $t$ and therefore the singularity initially moves
away from the real axis, regardless of the parameter values satisfying
$\epsilon\ll\alpha$. The high-order Fourier coefficients start increasing as
the singularity turns around, which occurs at $t\sim 0.37$, according to (71).


# Style modeling for Multi-Speaker Articulation-to-Speech

###### Abstract

In this paper, we propose a neural articulation-to-speech (ATS) framework that
synthesizes high-quality speech from articulatory signal in a multi-speaker
situation. Most conventional ATS approaches only focus on modeling contextual
information of speech from a single speaker’s articulatory features. To
explicitly represent each speaker’s speaking style as well as the contextual
information, our proposed model estimates style embeddings, guided from the
essential speech style attributes such as pitch and energy. We adopt
convolutional layers and transformer-based attention layers for our model to
fully utilize both local and global information of articulatory signals,
measured by electromagnetic articulography (EMA). Our model significantly
improves the quality of synthesized speech compared to the baseline in terms
of objective and subjective measurements in the Haskins dataset.

Index Terms—  Style-modeling, Multi-speaker, Speech synthesis, Articulation-
to-speech

##  1 Introduction

An articulation-to-speech (ATS) system aims to generate intelligible human
speech from articulatory movements. Recently, these systems have gained an
emerging popularity as an alternative communication tool for people having
difficulties in producing the intelligible speech [1] such as laryngectomees
[2]. As one of the essential parts of silent speech interface (SSI) [3], ATS
system can also bridge the connection between brain and speech signal [4].

Articulatory movements are represented as digital signals which are recorded
from places of articulation such as lips, tongue, and palate. Since these
kinematic recordings contain rich phonetic information of speech [5, 6], there
have been several statistical approaches that synthesizing human voice from
articulatory signals. Gaussian mixture model (GMM)-based ATS synthesizers are
proposed in [7, 8] to estimate Mel-Frequency Cepstral Coefficients (MFCCs)
from articulations. In [9], the effects of two different articulation styles
are investigated from synthetic acoustic feature from Hidden Markov Model
(HMM).

Thanks to recent advances in deep learning technologies, the performance of
ATS systems has significantly improved compared to the conventional signal
processing-based ones [10, 11, 12, 13, 14]. Nevertheless, there are some
limitations in those approaches. First, Most of the previous works synthesize
speech waveforms by solely utilizing the contextual information in
articulatory signals [1, 15]. Thus, there is still room for improvement of the
sound fidelity by providing a speaker’s representative acoustic features
(e.g., pitch and energy) to the synthesis process. Second, Most of the
previous ATS systems focus on a single-speaker synthesis scenario. Hence it
could be extended to a multi-speaker scenario with a unified system.

In this work, we propose a novel ATS system that estimates a target speaker’s
speaking style as well as contextual information from a given EMA signal. As
we explicitly model speaking styles of different speakers, our system can
support the multi-speaker voice generation scenario with a sole network. The
entire network architecture is composed of four modules: style encoder, mel
generator, and pitch and energy predictor. They are built on convolutional
neural networks (CNNs) and transformer-based attention layers to effectively
capture local and global characteristics from articulatory representations.
Specifically, the synthesis flow of the proposed model is as follows. First, a
style embedding is estimated from articulatory features with one-hot speaker
embedding, using a style encoder referenced by pitch and energy. The mel
generator then synthesizes mel-spectrum from EMA signal with the style
embedding. Finally, a neural vocoder generates raw waveform from the mel-
spectrum. By doing so, we can synthesize a highly intelligible speech with
speaker-specific characteristics from EMA signal under the multi-speaker
scenario. Audio samples are available online111https://miseulkim.github.io.

Our paper is organized as follows. In section 2, we summarize some previous
approaches related to our proposed model. We describe the detailed
architecture of our model in section 3. Experimental setup and results are
presented in section 4 and 5. Finally, we conclude the paper in section 6.

##  2 Related work

###  2.1 Neural articulation-to-speech (ATS) synthesis

Blessed by the development of deep neural networks, neural ATS systems have
significantly advanced performances. Those systems commonly estimate acoustic
features such as mel-spectra [1, 15], MFCCs [10], and vocoder parameters [13,
14] and use them as intermediate features for waveform synthesis. To construct
high-quality neural ATS systems, they encode high-dimensional latent
representations from input articulatory features using neural architectures
such as CNNs [15] or recurrent neural networks (RNNs) [1, 13, 14].

Introducing neural vocoders is also one of strategies to enhance the quality
of synthesized waveform [1, 10, 15]. These approaches achieve modeling of the
contextual information inherent in EMA signal, but there is still room for
more human-like speech. For this problem, [16] utilizes d-vector to represent
speaker information. However, the extracted speaker representation lacks time-
varying speaker characteristics because a time-averaging pooling is applied to
calculate the d-vector.

In our work, we also adopt the advanced neural architectures: CNNs,
Transformer, and neural vocoder to synthesize high-fidelity speech. Also, to
generate more natural speech, we estimate a speaker’s time-varying speaking
style from articulatory features and utilize it for the synthesis process.

###  2.2 Style modeling in text-to-speech (TTS)

Modeling of a speaker’s speaking style has been widely studied in TTS systems.
Fastspeech2 [17] utilizes speaker-specific prosody features such as pitch and
energy to improve the naturalness of synthesized speech. To explicitly model a
speaker’s speaking style as an embedding vector independently of text, [18]
proposes global style tokens (GSTs) and verifies their effectiveness on
various style-changing tasks including multi-speaker TTS. [19] uses a mel-
style encoder which is trained to capture speaker identity and prosody and
achieves a high-quality multi-speaker adaptive TTS system.

Similar to previous work, we explicitly model a speaker’s representative style
and apply our ATS system to the multi-speaker scenario. The difference from
aforementioned approaches is that we extract style information only from
articulatory features. We use a style encoder motivated from [19], which is
trained referencing additional acoustic features such as pitch and energy as
in [17].

##  3 Proposed Model

Fig. 1: Overall framework of the proposed model. For inference step, only the
part inside the light-yellow colored area is utilized.

In this section, we describe the architecture of the proposed model in detail.
As shown in Figure 1, our model consists of four modules: mel generator, style
encoder, pitch and energy predictor. The computation flow is given as follows.

First, to provide a target speaker’s information to perform a multi-speaker
ATS scenario, one-hot speaker embedding concatenated with EMA signal in
channel-wise is given as an input to the integration block. Next, this
combined vector passes through a CNN and gated-linear unit (GRU) layer with
residual connection (denoted as integration block). Then, the 256-dimensional
output vector of the integration block is fed into the style encoder to create
a style embedding. Continually, the obtained style embedding is used as an
input of the pitch and energy predictor in the training process and then pitch
and energy are extracted from. The style embedding is also concatenated with
the output embedding of the integration block to provide style information for
the mel generator to generate mel-spectrum. Finally, the generated mel-
spectrogram passes through the pre-trained neural vocoder to synthesize speech
waveform.

Table 1: Comparison of baseline and proposed model in MOS with confidence
intervals of 95%. Higher score indicated better quality

| Male speaker | Female speaker | Average
---|---|---|---
M01 | M02 | M03 | M04 | F01 | F02 | F03 | F04
Baseline | 2.82±0.11 | 2.21±0.11 | 2.81±0.12 | 2.75±0.12 | 2.70±0.11 | 2.58±0.11 | 2.75±0.11 | 2.71±0.11 | 2.67±0.04
Proposed | 3.25±0.12 | 2.63±0.11 | 3.27±0.12 | 3.33±0.11 | 3.35±0.11 | 3.12±0.12 | 3.62±0.11 | 3.36±0.11 | 3.24±0.04

###  3.1 Style encoder

EMA signals only contain a limited amount of information about the acoustic
properties of supralaryngeal configuratons (e.g., formants). However, human
voice is composed of various elements such as formant, pitch, energy and
prosody information depending on different characteristics of each speaker.
Therefore, we adopt a style encoder to obtain style vector which possesses an
additional acoustic information to generate expressive mel-spectrums from the
input EMA signals.

Based on the assumption that EMA is correlated with both contextual and
speaking style-related information, we directly estimates speech style
embedding from EMA by referencing pitch and energy component. Similar to the
structure of mel style encoder in [19], our style encoder consists of one
linear, two CNN blocks with residual connection, one multi-head attention and
one linear layer sequentially.

###  3.2 Pitch and Energy predictor

For pitch and energy predictor, we utilize the architecture of variance
predictor refers to [17], which consists of two CNN layers and one linear
layer. Each predictor generates $1*T$ frame-level pitch and energy output
vector respectively. To minimize the distance between the output of each
module and ground-truth target, L1 loss function is utilized for training.
Also, to reduce the training error, target pitch and energy are represented as
a decibel (dB) unit.

###  3.3 Mel generator

The mel generator is designed to generate high-quality mel-spectrum from a
given articulatory movement and a style embedding. For modeling both the
content-based global information and local features effectively from input
signal, we borrows the backbone of conformer [20], which integrates CNN and
transformer layer. For more details, our conformer block consists of one
convolution block with residual connection, transformer self-attention layer
with 4 head and two linear layers in order. Total 8 conformer blocks are
stacked for constructing generator.

In conclusion, total training loss function $L_{total}$ is defined as below
with weight parameters $\lambda_{mel}=0.8,\lambda_{pitch}=0.1$ and
$\lambda_{energy}=0.1$ respectively.

$L_{total}=\lambda_{mel}L_{mel}+\lambda_{pitch}L_{pitch}+\lambda_{energy}L_{energy}.$
(1)

As depicted in yellow-colored territory in Figure 1, only the integration
block, style encoder and mel-spectrum generator are used to generate mel-
spectrogram at the inference stage.

##  4 Experimental setup

###  4.1 Dataset

Haskins public dataset is utilized for our ATS task. The data consists of
correspondingly recorded EMA and audio file for 8 native American English
speakers. Average recording duration for per speaker is about 55 minutes.
Kinematic data and audio samples are recorded in 100 Hz and 44.1 kHz sampling
rate respectively. Among the recorded 8 point movements, 6 electrodes part are
used for our task, TR (tongue rear), TB (tongue blade), TT (tongue tip), UL
(upper lip), LL (lower lip), jaw (lower medial incisors). Each sensor position
includes trajectories orientation, X (posterior to anterior), Y (right to
left), Z (inferior to superior). Therefore, EMA vector of model’s input
consists of 18-dimensions (6 electrode * 3 dimension).

Table 2: STOI, PESQ, MCD and CER for each model. Summation symbol ’+’ means
that embedding components are added on the upper row system.

Model | STOI ↑ | PESQ ↑ | MCD ↓ | CER(%) ↓
---|---|---|---|---
Baseline | 0.663 | 1.154 | 6.809 | 31.14
Mel generator | 0.689 | 1.181 | 6.513 | 26.62
\+ Speaker embedding | 0.693 | 1.186 | 6.474 | 25.68
++ Style embedding | 0.697 | 1.193 | 6.436 | 26.01

###  4.2 Implementation details

Pre-processing Because kinematic data and audio samples are recorded in
substantially different sampling rates, we apply interpolation to EMA samples
for matching the number of sample as the one with target acoustic feature.
40-dimentional mel-spectrogram is the target for our approaches, which is
extracted with 1024 sample window length overlapping by 25%, considering the
low sampling rate of EMA. And the target audio is contaminated by background
noise sound, hence we remove the noise using speech enhancement module [21].
Also, we trim silence part from both of EMA and audio samples to prevent the
model confusion during training.

Baseline To perform a comparsion with SOTA transformer-based model, we build a
baseline [10], which stacks three residual convolution blocks and six
Transformer layers in order. Different from [10] which estimates MFCC
features, we train the baseline to estimate mel-spectrum from the articulatory
features. For training criteria, we adopt L1 loss to calculate the distance of
mel-spectrum between generated and ground truth.

Vocoder setup For obtaining the raw waveform from mel-spectrogram, hifi-GAN
[11] neural vocoder is adopted. Since the number of audio sample data in
haskin is insufficient for training, we first pre-train the hifi-GAN with
LJspeech dataset [22] and then fine-tune the model with haskins.

Fig. 2: Mel spectrum of reference and estimated samples for three systems.
(a): Target spectrum (b): Baseline (c): Mel generator with speaker embedding
(d): Mel generator with speaker embedding and style embedding.

###  4.3 Evaluation metric

We evaluate our results in terms of subjective and objective metrics.
Specifically, for subjective quality evaluation, we conduct MOS (mean opinion
score) test. For objective quality evaluation, we utilize four criteria, mel-
cepstral distortion (MCD), perceptual evaluation of speech quality (PESQ),
short-time objective intelligibility (STOI) and character error rate (CER). We
use a pre-trained ASR model— Wav2Vec [23] to compute CER.

##  5 Experimental results

###  5.1 Subjective quality measurements

To demonstrate the fidelity of the synthetic speech, we conduct a mean opinion
score (MOS) test between the baseline and proposed model. Total 20 respondents
are asked to give discrete scores from 1 to 5. We randomly select 20 sentences
from the test dataset per speaker and system.

As shown in the Table 1, our proposed network achieves much superb score than
the baseline for all speakers cross different genders. The experimental
results fully demonstrate the efficiency of proposed model of generating more
high-fidelity speech from articulatory signals than the baseline.

###  5.2 Objective quality measurements

We summarize objective quality of synthetic speech in Table 2. Our conformer-
based architectures show much superb performance than baseline model in terms
of MCD, PESQ, STOI. Experimental result shows that guiding the model with one-
hot speaker embedding could further improve the quality of the synthetic
speech. We also observe that adding both speaker and style embedding to the
network achieves the best performance among the networks. Thus, we conclude
that the two embeddings are definitely effective for enabling the network to
generate high-quality speech. However, we find that CER performance slightly
decline when a style embedding given to the network. It may result from the
synthetic speech conditioned on the style embedding has more prosody variation
compared to the baseline, therefore, it could suffer from the distortion in
phoneme characteristics [24].

Fig. 3: Mel spectrum, estimated pitch and energy trajectory for sentence ”The
wrist was badly strained and hung limp.”. (a): Target spectrum. (b): Synthetic
spectrum from proposed model. (c) Pitch trajectory. (d) Energy trajectory. The
blue line: target. Red line: estimated trajectory

###  5.3 Analysis of style embedding

We plot ground-truth and synthesized mel-spectrum obtained from three systems
in Figure 2 to confirm the quality of generated mel-sectrum. Spectrum in (b)
shows shorter length of formant compared with target spectrum. Spectrum in (c)
shows blurry boundary between formant. Spectrum (d) displays visually the
closest figure among three the ATS systems. To investigate the effectiveness
of style embedding more specifically, we visualize the estimated pitch and
energy in Figure 3. We can confirm that both the estimated pitch and energy
follow up reference well. From the results, we prove that the style embedding
exactly captures the speaker style information from the articulatory features.

##  6 Conclusion

In this paper, we introduced a novel high-quality speech synthesizer using
articulation signal (EMA) in a multi-speaker scenario. By estimating the style
embedding directly from EMA, which is definitely modeling pitch and spectrum
energy components, we generated not only intelligible but also speaker style
representative sound. Further fidelity of synthetic sound became increase by
informing the speaker information to the network using a speaker embedding. As
a model architecture, we adopted CNN and transformer attention layer based
style encoder and mel-spectrum generator to synthesize high-quality speech
from articulations. All experimental results proved that our proposed model
achieves higher performance in terms of both of objective and subjective
measurements.

##  7 Acknowledgements

This work was supported by the project ‘Alchemist Brain to X (B2X) Project’
through the Ministry of Trade, Industry and Energy (MOTIE), South Korea, under
Grant 20012355 and NTIS 1415181023.

## References

  * [1] Y.-W. Chen, K.-H. Hung, S.-Y. Chuang, J. Sherman, W.-C. Huang, X. Lu, and Y. Tsao, “Ema2s: An end-to-end multimodal articulatory-to-speech system,” in _2021 IEEE International Symposium on Circuits and Systems (ISCAS)_. IEEE, 2021, pp. 1–5.
  * [2] G. A. Gates, W. Ryan, J. Cooper Jr, G. F. Lawlis, E. Cantu, E. Lauder, R. W. Welch, and E. Hearne, “Current status of laryngectomee rehabilitation: I. results of therapy,” _American journal of otolaryngology_ , vol. 3, no. 1, pp. 1–7, 1982.
  * [3] B. Denby, T. Schultz, K. Honda, T. Hueber, J. M. Gilbert, and J. S. Brumberg, “Silent speech interfaces,” _Speech Communication_ , vol. 52, no. 4, pp. 270–287, 2010.
  * [4] G. K. Anumanchipalli, J. Chartier, and E. F. Chang, “Speech synthesis from neural decoding of spoken sentences,” _Nature_ , vol. 568, no. 7753, pp. 493–498, 2019.
  * [5] S. A. Kauffman, “Articulation of parts explanation in biology and the rational search for them,” in _Topics in the Philosophy of Biology_. Springer, 1976, pp. 245–263.
  * [6] L. E. Volaitis and J. L. Miller, “Phonetic prototypes: Influence of place of articulation and speaking rate on the internal structure of voicing categories,” _The Journal of the Acoustical Society of America_ , vol. 92, no. 2, pp. 723–735, 1992.
  * [7] T. Toda, A. W. Black, and K. Tokuda, “Mapping from articulatory movements to vocal tract spectrum with gaussian mixture model for articulatory speech synthesis,” in _Fifth ISCA Workshop on Speech Synthesis_ , 2004.
  * [8] S. Aryal and R. Gutierrez-Osuna, “Articulatory inversion and synthesis: towards articulatory-based modification of speech,” in _2013 IEEE International Conference on Acoustics, Speech and Signal Processing_. IEEE, 2013, pp. 7952–7956.
  * [9] B. Picart, T. Drugman, and T. Dutoit, “Continuous control of the degree of articulation in hmm-based speech synthesis,” in _Twelfth Annual Conference of the International Speech Communication Association_ , 2011.
  * [10] D. Gaddy and D. Klein, “An improved model for voicing silent speech,” in _Proceedings of the 59th Annual Meeting of the Association for Computational Linguistics and the 11th International Joint Conference on Natural Language Processing (Volume 2: Short Papers)_. Online: Association for Computational Linguistics, Aug. 2021, pp. 175–181. [Online]. Available: https://aclanthology.org/2021.acl-short.23
  * [11] J. Kong, J. Kim, and J. Bae, “Hifi-gan: Generative adversarial networks for efficient and high fidelity speech synthesis,” _Advances in Neural Information Processing Systems_ , vol. 33, pp. 17 022–17 033, 2020.
  * [12] B. Cao, A. Wisler, and J. Wang, “Speaker adaptation on articulation and acoustics for articulation-to-speech synthesis,” _Sensors_ , vol. 22, no. 16, p. 6056, 2022.
  * [13] B. Cao, M. J. Kim, J. R. Wang, J. P. van Santen, T. Mau, and J. Wang, “Articulation-to-speech synthesis using articulatory flesh point sensors’ orientation information.” in _INTERSPEECH_ , 2018, pp. 3152–3156.
  * [14] F. Taguchi and T. Kaburagi, “Articulatory-to-speech conversion using bi-directional long short-term memory.” in _Interspeech_ , 2018, pp. 2499–2503.
  * [15] T. G. Csapó, C. Zainkó, L. Tóth, G. Gosztolya, and A. Markó, “Ultrasound-based articulatory-to-acoustic mapping with waveglow speech synthesis,” _arXiv preprint arXiv:2008.03152_ , 2020.
  * [16] K. Katsurada and K. Richmond, “Speaker-Independent Mel-Cepstrum Estimation from Articulator Movements Using D-Vector Input,” in _Proc. Interspeech 2020_ , 2020, pp. 3176–3180.
  * [17] Y. Ren, C. Hu, X. Tan, T. Qin, S. Zhao, Z. Zhao, and T.-Y. Liu, “Fastspeech 2: Fast and high-quality end-to-end text to speech,” _arXiv preprint arXiv:2006.04558_ , 2020.
  * [18] Y. Wang, D. Stanton, Y. Zhang, R.-S. Ryan, E. Battenberg, J. Shor, Y. Xiao, Y. Jia, F. Ren, and R. A. Saurous, “Style tokens: Unsupervised style modeling, control and transfer in end-to-end speech synthesis,” in _International Conference on Machine Learning_. PMLR, 2018, pp. 5180–5189.
  * [19] D. Min, D. B. Lee, E. Yang, and S. J. Hwang, “Meta-stylespeech: Multi-speaker adaptive text-to-speech generation,” in _International Conference on Machine Learning_. PMLR, 2021, pp. 7748–7759.
  * [20] A. Gulati, J. Qin, C.-C. Chiu, N. Parmar, Y. Zhang, J. Yu, W. Han, S. Wang, Z. Zhang, Y. Wu _et al._ , “Conformer: Convolution-augmented transformer for speech recognition,” _arXiv preprint arXiv:2005.08100_ , 2020.
  * [21] T. Sainburg, M. Thielk, and T. Q. Gentner, “Finding, visualizing, and quantifying latent structure across diverse animal vocal repertoires,” _PLoS computational biology_ , vol. 16, no. 10, p. e1008228, 2020.
  * [22] K. Ito and L. Johnson, “The lj speech dataset,” https://keithito.com/LJ-Speech-Dataset/, 2017.
  * [23] A. Baevski, Y. Zhou, A. Mohamed, and M. Auli, “wav2vec 2.0: A framework for self-supervised learning of speech representations,” _Advances in Neural Information Processing Systems_ , vol. 33, pp. 12 449–12 460, 2020.
  * [24] N. U. Rani and P. Girija, “Error analysis to improve the speech recognition accuracy on telugu language,” _Sadhana_ , vol. 37, no. 6, pp. 747–761, 2012\.