HeunD (Double Confluent Heun Function)

Note

This page focuses on the double confluent Heun equation (DCHE) and the standard local confluent Heun function commonly called HeunD in computer algebra systems.

What HeunD is

Heun’s general equation (four regular singular points) admits several confluent limits in which singularities merge. The double confluent limit is the most singular of the classical confluences: two pairs of regular singularities coalesce so that the resulting equation has

• exactly two singular points (at $z=0$ and $z=\infty$),
• and both are irregular (rank $1$ in the generic case).

This puts DCHE in the same conceptual family as equations studied with Stokes phenomena, connection coefficients, and Floquet theory. In physics, DCHE frequently appears after variable changes that convert a radial or angular ODE into something with exponential behavior at both ends, and it also contains classical periodic problems (e.g. Mathieu-type equations) as special cases.

---

Conventions used on this page

Always write the ODE and the normalization

“HeunD” is not a single universal object: different sources use different canonical equations, different parameter names, and even different normalization points.
To be unambiguous, always state: (i) the exact ODE you mean and (ii) the normalization (e.g. $y(1)=1$, $y'(1)=0$).

You will see (at least) four conventions in real research workflows:

1. DLMF canonical DCHE (NIST Digital Library of Mathematical Functions).
2. Wolfram / Mathematica HeunD convention.
3. Maple HeunD convention (different singularity locations).
4. Ronveaux (1995) / “symmetric canonical form” used in Heun’s Differential Equations (Part C).

This page uses Wolfram/DLMF for the computational definition of HeunD, and Ronveaux (1995) for asymptotics and connection intuition, because Part C gives a clean asymptotic basis at both irregular singularities.

To keep symbols from colliding, we’ll label parameter sets explicitly:

• Wolfram parameters: $(q,\alpha,\gamma,\delta,\epsilon)$.
• DLMF parameters: $(q,\alpha,\gamma_{\rm D},\delta_{\rm D})$ (DLMF fixes one scale).
• Ronveaux parameters: $(\alpha_R,\beta_{-1},\beta_1,\gamma_R)$ plus the operator $D=z\,d/dz$.
• Maple parameters: $(\alpha_M,\beta_M,\gamma_M,\delta_M)$.

---

1. Canonical forms of the double confluent Heun equation

1.1 DLMF standard form (two rank-1 irregular singularities)

A standard canonical DCHE used by DLMF is

$$y''(z) +\left(\frac{\delta_{\rm D}}{z^2}+\frac{\gamma_{\rm D}}{z}+1\right) y'(z) +\frac{\alpha z-q}{z^2}\,y(z)=0.$$

• Singularities: irregular at $z=0$ and $z=\infty$, each of rank 1.
• Any $z_0\neq 0$ is an ordinary point, so local solutions are analytic around $z_0$.

1.2 Wolfram/Mathematica form and the definition of HeunD

Wolfram’s documentation uses the equivalent polynomial-coefficient form

$$z^2 y''(z) + (\gamma+\delta z+\epsilon z^2) y'(z) + (\alpha z-q) y(z)=0,$$

or, dividing by $z^2$,

$$y''(z)+\left(\frac{\gamma}{z^2}+\frac{\delta}{z}+\epsilon\right) y'(z)+\frac{\alpha z-q}{z^2}y(z)=0.$$

Normalization (Wolfram’s HeunD). The Wolfram Language defines HeunD[q,α,γ,δ,ϵ,z] as the power-series solution of the above ODE satisfying

$$y(1)=1,\qquad y'(1)=0.$$

So, in Wolfram notation,

$$\mathrm{HeunD}(q,\alpha,\gamma,\delta,\epsilon;z) \equiv y(z)\ \text{with}\ y(1)=1,\ y'(1)=0.$$

Why $z=1$?

In this canonical form, $z=0$ is an irregular singular point, so a Frobenius normalization at $z=0$ is not available.
Choosing $z_0=1$ fixes a unique analytic branch in the disk $|z-1|<1$ and then extends it by analytic continuation.

Relation to DLMF.

• If $\epsilon\neq 0$, the rescaling $z=t/\epsilon$ reduces the equation to a form with $\epsilon=1$. Concretely, with $y(z)=Y(t)$ and $t=\epsilon z$, you get

  $$t^2 Y''(t) + (\epsilon\gamma+\delta t+t^2)Y'(t) + \left(\frac{\alpha}{\epsilon}t-q\right)Y(t)=0,$$

  so (up to notation) $\epsilon$ can be set to 1 by the parameter redefinitions

  $$\gamma\mapsto \epsilon\gamma,\qquad \alpha\mapsto \alpha/\epsilon,\qquad \delta\mapsto \delta,\qquad q\mapsto q.$$

• After setting $\epsilon=1$, DLMF and Wolfram differ only by a swap of parameter names in the $1/z^2$ and $1/z$ coefficients:

$$\gamma_{\rm D} = \delta,\qquad \delta_{\rm D}=\gamma,\qquad (\epsilon=1).$$

1.3 Ronveaux (1995) symmetric canonical form (Part C)

Part C of Heun’s Differential Equations works with the “symmetric canonical” DCHE

$$D^2 y +\alpha_R\left(z+\frac{1}{z}\right) D y +\left( \left(\beta_1+\frac12\right)\alpha_R z +\left(\frac{\alpha_R^2}{2}-\gamma_R\right) +\left(\beta_{-1}-\frac12\right)\frac{\alpha_R}{z} \right)y =0,$$

where

$$D := z\frac{d}{dz}.$$

In this form, $z=0$ and $z=\infty$ are again rank-1 irregular singularities, and the coefficients are meromorphic on $\mathbb{C}^*=\mathbb{C}\setminus\{0\}$.

A key structural fact used repeatedly in Part C is that the DCHE is stable under gauge transformations of the form

$$y(z)=z^{\zeta_0}\exp(\zeta_1 z-\zeta_{-1}/z)\,u(z),$$

together with simple changes of variable ($z\mapsto \tau z$ and $z\mapsto \tau/z$). These are the natural operations for an equation whose singularities are irregular.

1.4 Maple’s HeunD convention (irregular singularities at $\pm 1$)

Maple’s built-in HeunD(α,β,γ,δ,x) solves a different but equivalent normal form of the DCHE, with irregular singularities placed symmetrically at $x=\pm 1$ and the origin $x=0$ as a regular point.

One convenient explicit statement of Maple’s equation is (see Appendix A.5 in Birkandan 2020)

$$y''(x) - \frac{\alpha_M x^4 - 2 x^5 + 4 x^3 - \alpha_M - 2 x}{(x-1)^3(x+1)^3}\,y'(x) - \frac{-x^2\beta_M+(-\gamma_M-2\alpha_M)x-\delta_M}{(x-1)^3(x+1)^3}\,y(x)=0,$$

and Maple computes HeunD as a power series around $x=0$ (radius of convergence 1, limited by $x=\pm 1$).

Do not mix Maple and Wolfram parameters

Even when both systems call the function HeunD, the parameters and even the location of singularities differ.
When comparing results across CAS, always compare by substituting back into the ODE.

---

2. What “the function” HeunD really is

Because the DCHE has no regular singular points, the clean “Frobenius at $z=0$” definition that works for HeunG and HeunC is not available.

Mathematically, a “double confluent Heun function” is best thought of as:

• choose a canonical DCHE and a base point $z_0$ that is an ordinary point (e.g. $z_0=1$ in Wolfram’s convention),
• define the unique local analytic solution by prescribing two initial conditions (normalization),
• extend by analytic continuation along a path that avoids the singular point(s).

This automatically implies two practical facts:

1. HeunD is branch dependent. Looping around $z=0$ changes the solution by monodromy (Part C discusses this in detail).
2. You cannot assign a single “value at $z=0$” in general. In the Wolfram/DLMF form, $z=0$ is an irregular singular point, and the value is not well-defined as a regular limit.

---

3. Local power series: the defining expansion used in CAS

3.1 Wolfram HeunD: series around $z=1$

Let $t=z-1$ and expand

$$y(z)=\sum_{n=0}^\infty c_n\,t^n, \qquad c_0=y(1)=1,\quad c_1=y'(1)=0.$$

Plugging this into Wolfram’s canonical ODE

$$z^2 y'' + (\gamma+\delta z+\epsilon z^2)y' + (\alpha z-q)y=0$$

gives a linear recurrence for the Taylor coefficients. Define

$$A_0:=\gamma+\delta+\epsilon,\qquad A_1:=\delta+2\epsilon,\qquad A_2:=\epsilon,\qquad B_0:=\alpha-q,\qquad B_1:=\alpha.$$

Then the coefficients satisfy, for $n\ge 0$ (with the convention $c_{-1}=0$),

$$(n+2)(n+1)c_{n+2} +(n+1)(2n+A_0)c_{n+1} +\big(n(n-1)+A_1 n + B_0\big)c_n +\big(A_2(n-1)+B_1\big)c_{n-1} =0.$$

From $c_0=1$, $c_1=0$ one finds the first nontrivial terms

$$c_2=\frac{q-\alpha}{2},\qquad c_3=-\frac{(2+A_0)(q-\alpha)+\alpha}{6},$$

and so on.

Convergence radius. In this canonical form the nearest singularity to $z=1$ is $z=0$, so the Taylor series converges at least for $|z-1|<1$. Outside this disk you must use analytic continuation (CAS does this internally) or solve the ODE numerically along a path.

3.2 Maple HeunD: series around $x=0$

In Maple’s convention the ordinary point is $x=0$ (singularities at $x=\pm 1$), and HeunD(α,β,γ,δ,x) is computed from the Taylor series around $x=0$ with radius of convergence 1.

If you want to reproduce Maple’s series yourself, the workflow is identical:

1. write $y(x)=\sum c_n x^n$,
2. plug into the Maple DCHE,
3. solve the resulting recurrence for $c_n$ using the chosen normalization.

---

4. Asymptotic solutions and Stokes phenomenon (Ronveaux form)

When your application imposes boundary conditions near $z=0$ and/or $z=\infty$, you need asymptotics, not just Taylor series.

Part C constructs canonical asymptotic bases at both irregular singularities for the symmetric canonical DCHE (Section 2.2).

4.1 As $z\to\infty$: two canonical behaviours

There exist two independent solutions $y_{\infty 1}(z)$ and $y_{\infty 2}(z)$ with characteristic behaviours

$$y_{\infty 1}(z) \sim (\alpha_R z)^{-\beta_1-\tfrac12} \left(1+O\left(\frac{1}{z}\right)\right),$$

and

$$y_{\infty 2}(z) \sim e^{-\alpha_R z}(\alpha_R z)^{\beta_1-\tfrac12} \left(1+O\left(\frac{1}{z}\right)\right),$$

valid in appropriate Stokes sectors (the precise angular ranges are given in Theorem 2.2.14 of Part C).

So one solution is essentially power-like, while the other is exponentially small (or large, depending on $\Re(\alpha_R z)$) at infinity.

4.2 As $z\to 0$: two canonical behaviours

Similarly, there exist $y_{01}(z)$ and $y_{02}(z)$ with

$$y_{01}(z) \sim e^{\alpha_R/z}\left(\frac{\alpha_R}{z}\right)^{-\beta_{-1}-\tfrac12} \left(1+O(z)\right),$$

and

$$y_{02}(z) \sim \left(\frac{\alpha_R}{z}\right)^{\beta_{-1}-\tfrac12} \left(1+O(z)\right),$$

again in appropriate sectors (Theorem 2.2.20 of Part C).

4.3 Asymptotic coefficient recurrences (useful in practice)

Part C also gives a practical way to compute the coefficients of the asymptotic series at infinity. Writing

$$y_{\infty 1}(z) \sim (\alpha_R z)^{-\beta_1-\tfrac12} \sum_{n=0}^\infty \psi_n (\alpha_R z)^{-n}, \qquad \psi_0=1,$$

the coefficients satisfy a three-term recurrence (Eq. 2.2.6 in Part C):

$$(n+1)\psi_{n+1} -\left(\left(n+\beta_1+\tfrac12\right)^2-\gamma_R+\frac{\alpha_R^2}{2}\right)\psi_n +\alpha_R^2\left(n+\beta_1-\beta_{-1}\right)\psi_{n-1} =0, \qquad n\in\mathbb{N},$$

with $\psi_{-1}=0$.

How this helps computations

If your numerical evaluation is unstable near $z=\infty$ (or $z=0$), use the corresponding asymptotic series to generate boundary data in a safe sector, then integrate the ODE toward the region of interest.

4.4 A Wronskian identity you can use as a check

For the symmetric canonical DCHE, it is convenient to use the $D$-Wronskian

$$W_D[y_1,y_2](z):= y_1(z)\,D y_2(z) - y_2(z)\,D y_1(z),\qquad D=z\frac{d}{dz}.$$

Part C shows that for any two solutions $y_1,y_2$ of the symmetric canonical equation,

$$W_D[y_1,y_2](z)=\mathrm{const}\cdot \exp\left(-\alpha_R z+\frac{\alpha_R}{z}\right).$$

This is a practical numerical sanity check: if you integrate the DCHE numerically and compute $W_D[y_1,y_2](z)\,\exp(\alpha_R z-\alpha_R/z)$ along your path, the result should be (approximately) constant.

---

5. Connection problems and the “accessory parameter” as an eigenvalue

In many physics applications, the DCHE does not appear with all parameters fixed. Instead:

• some parameters are fixed by local exponents or physical constants,
• one parameter (often called an accessory parameter, e.g. $q$ in DLMF/Wolfram or $\gamma_R$ in Ronveaux’s form) is determined by imposing boundary conditions at both irregular singularities.

Part C formulates this via a connection matrix $C$ that relates the asymptotic bases at $0$ and $\infty$:

$$Y_0(z) = Y_{\infty}(z)\,C,$$

where $Y_0=(y_{01},y_{02})$ and $Y_{\infty}=(y_{\infty 1},y_{\infty 2})$ are fundamental solution matrices built from the canonical asymptotic solutions.

A boundary condition like “select the subdominant solution at $z\to 0$ and at $z\to\infty$” is then a condition that a specific entry of $C$ vanishes (equivalently: a Wronskian/connection coefficient is zero).

5.1 A practical numerical workflow for eigenvalues

If your goal is to determine the accessory parameter (say $q$):

1. Choose a canonical form and map your ODE to it.
2. Choose sectors where the desired asymptotics are stable (subdominant).
3. For a trial $q$, evaluate the normalized local solution (e.g. the Wolfram HeunD) and match it to asymptotic data near $0$ and/or $\infty$.
4. Compute the unwanted mixing coefficient (connection coefficient).
5. Use a root finder in $q$.

5.2 A concrete special value: $q$ at $\alpha=0$

A particularly instructive checkpoint is the degenerate limit $\alpha=0$. In that limit the connection coefficient has a closed form in terms of Gamma functions:

$$q(0,\beta_{-1},\beta_1,\gamma) = -2\pi e^{-i\pi\beta_1}\, \frac{1}{ \Gamma\!\left(\frac12+\beta_1+\gamma\right) \Gamma\!\left(\frac12+\beta_1-\gamma\right)}.$$

---

6. Symmetries and transformations you can exploit

Because DCHE coefficients are rational (or meromorphic) and the singularities are irregular, gauge and variable transformations are often the fastest way to make computations stable.

Two families are particularly important (see Part C, Section 1.1–1.2):

1. Gauge transforms

   $$y(z)=z^{\zeta_0}e^{\zeta_1 z-\zeta_{-1}/z}\,u(z),$$

   which can simplify the dominant exponential behaviours at $0$ and $\infty$.

2. Inversion/scaling

   $$z\mapsto \tau z,\qquad z\mapsto \tau/z,$$

   which swap or rescale the singular points.

These transformations underlie most “identity” formulas you see in CAS and in the literature, and they are also a practical tool: map the point where you need a value to a region where the Taylor series converges rapidly.

---

7. Degenerate and special cases (reductions)

The DCHE contains many classical equations as special/degenerate limits. Part C (Section 1.3) gives a classification in terms of invariants of the general DCHE, with reductions to

• Mathieu-type equations (after $z=e^{ix}$ in symmetric normal-form settings),
• Bessel equations (when only one side of the Laurent structure survives),
• confluent hypergeometric-type equations,
• and Euler-type equations

in appropriate parameter regimes.

From a computational perspective, these reductions matter because:

• they provide sanity checks (compare HeunD against Bessel/Mathieu in the limit),
• they often give better-conditioned representations for extreme parameter values.

---

8. Floquet solutions and convergent Laurent expansions on $\mathbb C^*$

A defining feature of DCHE (and one of the most useful computational handles) is the existence of solutions with controlled monodromy around $0$.

8.1 Monodromy operator and characteristic exponent

Let $m$ denote analytic continuation around $z=0$ once:

$$(my)(z)=y(e^{2\pi i}z).$$

A Floquet solution is a solution satisfying

$$y(e^{2\pi i}z)=e^{2\pi i\nu}\,y(z),$$

for some $\nu\in\mathbb C$ called the characteristic exponent (or Floquet exponent).

Because we are working on $\Omega$ (the log surface), the expression $z^\nu=e^{\nu\log z}$ is well-defined there.

8.2 Laurent-series representation

Ronveaux shows that a Floquet solution can be represented as

$$y(z)=z^\nu h(z),$$

where $h$ is holomorphic on $\mathbb C^*$ and therefore has a convergent Laurent series

$$h(z)=\sum_{n=-\infty}^{\infty}\eta_n z^n \qquad (z\in\mathbb C^*).$$

Substituting $y=z^\nu\sum_{n\in\mathbb Z}\eta_n z^n$ into the symmetric DCHE yields a three-term recurrence valid for all $n\in\mathbb Z$:

$$\alpha\left(n+\nu+\beta_{-1}+\frac12\right)\eta_{n+1} +\left((n+\nu)^2-\gamma+\frac{\alpha^2}{2}\right)\eta_n +\alpha\left(n+\nu+\beta_1-\frac12\right)\eta_{n-1} =0, \qquad n\in\mathbb Z.$$

This recurrence is central for:

• computing Floquet solutions numerically via continued fractions / minimal-solution conditions,
• studying when $\nu$ is integral or half-integral (single-valuedness or sign changes),
• and formulating eigenvalue problems in $\gamma$.

Caution

The recurrence is two-sided ($n\in\mathbb Z$). Picking a solution that converges/behaves correctly in both $n\to+\infty$ and $n\to-\infty$ is a nontrivial constraint and is exactly where Floquet/continued-fraction methods enter.

---

9. Quasi-polynomials: when the “asymptotics” actually terminate

A striking phenomenon in DCHE is the existence of special solutions whose asymptotic series terminate in a way that produces a polynomial factor—analogous in spirit to Heun polynomials in the regular case.

Ronveaux calls these quasi-polynomials. A quasi-polynomial solution of degree $N$ has the form

$$y(z)=z^\nu \exp\!\left(-\alpha\left(\varepsilon_+ z+\varepsilon_-\frac1z\right)\right) P_N(z), \qquad z\in\Omega,$$

where

• $\nu\in\mathbb C$,
• $\varepsilon_\pm\in\{0,1\}$,
• and $P_N$ is a polynomial of degree $N$ with $P_N(0)\neq0$.

A particularly important subcase is the “pure polynomial” Floquet factor:

$$y(z)=z^\nu P_N(z).$$

Ronveaux shows that such a solution exists only under discrete constraints linking $(\beta_{-1},\beta_1)$ and a termination condition on the asymptotic coefficients $\psi_n$ of the canonical $\Psi$:

• one must have $\beta_{-1}-\beta_1=N+1$ (with $N\in\mathbb N^*$),
• and the coefficient $\psi_{N+1}$ must vanish.

Since $\psi_{N+1}$ is (for fixed $\alpha,\beta_{-1},\beta_1$) a polynomial in $\gamma$ of degree $N+1$, the condition $\psi_{N+1}=0$ yields a finite set of candidate accessory parameters $\gamma$ that produce quasi-polynomials.

Tip

Quasi-polynomials are the DCHE analogue of “exactly solvable” parameter points: they convert a global analytic-connection problem into an algebraic condition in $\gamma$.

---

10. DCHE as a Floquet eigenvalue problem

A powerful point of view—especially for physics—is to treat $\gamma$ (or a linear shift of it) as an eigenvalue, with the Floquet exponent $\nu$ playing the role of the quasi-momentum/Bloch phase.

Ronveaux formulates an eigenvalue problem using the operator

$$G(\beta)y := z\left(D+\beta_1+\frac12\right)y+\frac1z\left(D+\beta_{-1}-\frac12\right)y,$$

and considers

$$\left(D^2+\alpha\,G(\beta)-\lambda\right)y=0,$$

with the identification

$$\lambda=-\frac{\alpha^2}{2}+\gamma.$$

For fixed $(\nu,\alpha,\beta_{-1},\beta_1)$, the eigenvalues form a countable family $\lambda_\mu(\alpha,\beta)$ indexed by $\mu\in\nu+\mathbb Z$, and in the unperturbed case $\alpha=0$ one has

$$\lambda_\mu(0,\beta)=\mu^2.$$

Moreover, for small $|\alpha|$ (in an appropriate domain), $\lambda_\mu(\alpha,\beta)$ is holomorphic in $\alpha^2$ and admits a power series expansion

$$\lambda_\mu(\alpha,\beta)=\mu^2+\sum_{m=1}^\infty \lambda_{\mu m}(\beta)\,\alpha^{2m}.$$

The first perturbative coefficient is explicitly

$$\lambda_{\mu 1}(\beta)= -\frac12+\frac{2\beta_{-1}\beta_1}{4\mu^2-1}.$$

Note

This is one of the cleanest “graduate-student entry points” into DCHE: start from $\alpha=0$ where everything is elementary, then view DCHE as a deformation controlled by $\alpha$, with eigenvalues expanding in $\alpha^2$.

---

References and further reading

• A. Ronveaux (ed.), Heun’s Differential Equations, Oxford University Press (1995). Part C: Double Confluent Heun Equation.

• NIST Digital Library of Mathematical Functions (DLMF), §31.12 “Confluent Forms of Heun’s Equation” (see Eq. 31.12.2 for the doubly confluent equation): https://dlmf.nist.gov/31.12. The equation image for 31.12.2 is at https://dlmf.nist.gov/31.12.E2.

• Maple Help: HeunD (double confluent Heun function): https://www.maplesoft.com/support/help/maple/view.aspx?path=HeunD.

• Wolfram Language documentation: HeunD (https://reference.wolfram.com/language/ref/HeunD.html) and HeunDPrime (https://reference.wolfram.com/language/ref/HeunDPrime.html).

• T. Birkandan, “An algorithm for the analytic solution of second order linear differential equations” (2020), Appendix A: https://arxiv.org/abs/2010.01563.