# Shift of Dirac points and strain induced pseudo-magnetic field in graphene

Hua Tong Yang<EMAIL_ADDRESS>Center for Advanced Optoelectronic
Functional Materials Research, Key Laboratory for UV-Emitting Materials and
Technology of Ministry of Education, and School of Physics, Northeast Normal
University, Changchun 130024, China

###### Abstract

We propose that the strain induced effective pseudo-magnetic field in graphene
can also be explained by a curl movement of the Dirac points, if the Dirac
points can be regarded as a slowly varying function of position. We also prove
that the Dirac points must be confined within two triangles, each one has 1/8
the area of the Brillouin zone.

###### pacs:

73.22.Pr, 73.22.Dj, 73.22.Gk, 73.20.At

The discovery of graphene, a monolayer carbon atom sheet Novoselov-sci306 ,
and the development of experimental technique to manipulate this two-
dimensional(2D) material have ignited intense interest in this system Geim ;
Meyer ; Vozmediano ; Neto . One of the most attractive characters of graphene
is that its low energy excitation satisfies a massless 2D Dirac equation
Semenoff , and the chemical potential crosses its Dirac points(or Fermi
points) in neutral graphene. These special characters lead to many unusual
properties and new phenomena Gusynin ; Neto ; Novoselov2005 ; Zhang2005 , such
as the anomalous integer quantum Hall effect(QHE) Novoselov2005 ; Zhang2005 .
Recently, experiments have confirmed another remarkable effect that mechanical
strain can induce a very strong effective pseudo-magnetic field, leading to a
pseudo-QHE, which can be observed in zero magnetic field Guinea-10 ; Levy . In
this paper we propose that the strain induced effective vector potential can
also be explained by shift $\delta\mathbf{K}(\mathbf{x})$ of the Dirac points
$\mathbf{K}(\mathbf{x})$, its effective pseudo-magnetic field is in proportion
to $\nabla\times\mathbf{K}(\mathbf{x})$, only if the Dirac points
$\mathbf{K}(\mathbf{x})$ can be regarded as a slowly varying function of
position, and the Fermi velocity is generalized to a tensorZhu-Wang . We also
prove that the Dirac points can not be arbitrarily moved, they must be
confined within two triangles, each one has 1/8 the area of the Brillouin
zone(BZ).

Firstly, consider a tight-binding Hamiltonian describing a uniformly deformed
honeycomb lattice with three different nearest-neighbor hopping energies
$t_{1},t_{2},t_{3}$Hasegawa ; Montambaux08 ; Montambaux09 :

$\displaystyle\hat{H}=-\sum_{<\mathbf{i}a,\mathbf{j}b>}t_{\mathbf{i}a,\mathbf{j}b}c^{{\dagger}}_{\mathbf{i}a}c_{\mathbf{j}b}+h.c.,$
(1)

where $c_{\mathbf{j}b}$ ($c^{{\dagger}}_{\mathbf{i}a}$) are
annihilation(creation) operators, $\mathbf{i}$($\mathbf{j}$) are position
vectors of unit cells, $a$($b$) denote two inequivalent atoms in a unit cell,
$t_{\mathbf{i}a,\mathbf{j}b}$ is the electronic hopping energy from the
$\mathbf{j}$th unit cell $b$ atom to $\mathbf{i}$th unit cell $a$ atom.
Suppose that the deformed lattice remains invariant under spatial translation,
i.e., $t_{\mathbf{i}a,\mathbf{j}b}$ only depends on ${\bf i}-{\bf j}$, but the
three nearest-neighbor hopping energies $t_{1,2,3}$ may be different owing to
anisotropy of strains, as shown in Fig.1.

Figure 1: Unit cell and hopping parameters for deformed graphene.

The hopping parameters can be written as some $2\times 2$ matrixes
$\mathbf{t}(\mathbf{i}-\mathbf{j})$, whose elements are defined by $[{\bf
t}(\mathbf{i}-\mathbf{j})]_{a,b}\equiv t_{\mathbf{i}a,\mathbf{j}b}.$ For this
nearest-neighbor tight-binding Hamiltonian, the non-vanishing hopping matrixes
are

$\displaystyle{\bf t}(0)=\bigg{(}\begin{array}[]{cc}0&t_{1}\\\
t_{1}&0\end{array}\bigg{)},{\bf t}({\bf
a}_{1})=\bigg{(}\begin{array}[]{cc}0&t_{2}\\\ 0&0\end{array}\bigg{)},{\bf
t}({\bf a}_{2})=\bigg{(}\begin{array}[]{cc}0&t_{3}\\\ 0&0\end{array}\bigg{)}$
(8)

and ${\bf t}(-{\bf a}_{1})={\bf t}^{\dagger}({\bf a}_{1}),{\bf t}(-{\bf
a}_{2})={\bf t}^{\dagger}({\bf a}_{2}).$ By Fourier transformation

$c_{\mathbf{j},a(b)}=\frac{1}{\sqrt{N}}\sum_{\mathbf{k}}c_{\mathbf{k},a(b)}\exp(i\mathbf{k}\cdot\mathbf{j})$

with $N$ a normalization constant, the Hamiltonian (1) can be cast into the
form

$\displaystyle\hat{H}=-\sum_{\mathbf{k}}\big{[}c^{{\dagger}}_{\mathbf{k},a},c^{{\dagger}}_{\mathbf{k},b}\big{]}\bigg{[}\begin{array}[]{cc}h_{aa}(\mathbf{k})&h_{ab}(\mathbf{k})\\\
h_{ba}(\mathbf{k})&h_{bb}(\mathbf{k})\end{array}\bigg{]}\bigg{[}\begin{array}[]{c}c_{\mathbf{k},a}\\\
c_{\mathbf{k},b}\end{array}\bigg{]},$ (13)

where $h_{aa}(\mathbf{k})=h_{bb}(\mathbf{k})=0,$
$h_{ab}(\mathbf{k})=h_{ba}^{\ast}(\mathbf{k}),$ and

$\displaystyle
h_{ba}(\mathbf{k})=t_{1}+t_{2}\exp(i\mathbf{k}\cdot\mathbf{a}_{1})+t_{3}\exp(i\mathbf{k}\cdot\mathbf{a}_{2}),$
(14)

with $\mathbf{a}_{1},$$\mathbf{a}_{2}$ the lattice unit vectors. The energy
bands obtained by diagonalizing this Hamiltonian areWallace

$E_{\pm}(\mathbf{k})=\pm|t_{1}+\tilde{t}_{2}(\mathbf{k})+\tilde{t}_{3}(\mathbf{k})|,$
(15)

where $\tilde{t}_{2}(\mathbf{k})=t_{2}e^{i\mathbf{k}\cdot\mathbf{a}_{1}},$
$\tilde{t}_{3}(\mathbf{k})=t_{3}e^{i\mathbf{k}\cdot\mathbf{a}_{2}}$, the plus
sign corresponds to the upper($\pi$) and minus to the lower($\pi^{*}$) band
respectively. From Eq.(15) we notice that if $\mathbf{K}$ is a zero point of
$h_{ba}(\mathbf{K})$, i.e.,

$t_{1}+\tilde{t}_{2}(\mathbf{K})+\tilde{t}_{3}(\mathbf{K})=0,$ (16)

then $E_{+}(\mathbf{k})$ and $E_{-}(\mathbf{k})$ will meet at $\mathbf{K}$,
i.e., $E_{+}(\mathbf{K})=E_{-}(\mathbf{K})=0$, this $\mathbf{K}$ is known as
the Dirac point. The Hamiltonian (13) can be expanded up to a linear order in
$\mathbf{p}=\mathbf{k}-\mathbf{K}$ in a neighborhood of point $\mathbf{K}$

$\displaystyle\bigg{[}\begin{array}[]{cc}h_{aa}(\mathbf{k})&h_{ab}(\mathbf{k})\\\
h_{ba}(\mathbf{k})&h_{bb}(\mathbf{k})\end{array}\bigg{]}\simeq\bigg{[}\begin{array}[]{cc}0&\vec{\alpha}^{*}\cdot\mathbf{p}\\\
\vec{\alpha}\cdot\mathbf{p}&0\end{array}\bigg{]}=v_{\mu\nu}\sigma^{\mu}p^{\nu},$
(21)

where $\mu,\nu=1,2$ denote two components of a 2D vector and a sum over the
repeated indices $\mu,\nu$ is implied, $\vec{\alpha}$ is a complex vector with
$\textrm{Re}(\vec{\alpha})=(v_{11},v_{12})$,
$\textrm{Im}(\vec{\alpha})=(v_{21},v_{22})$, $\sigma^{1,2}$ are Pauli matrixes
acting on the sublattice degree of freedom, tensor $v_{\mu\nu}$ represents the
anisotropy of the dispersion near the Dirac points, it only occurs noticeable
departure from $v_{F}\delta_{\mu\nu}$ in a strongly deformed grapheneZhu-Wang
. However, after this modification the strain induced effective vector
potential will acquire a direct physical meaning. For a graphene under
nonuniform but slowly varying strain, $t_{i}(\mathbf{x})$ and hence the Dirac
point $\mathbf{K}(\mathbf{x})$ as well as $v_{\mu\nu}(\mathbf{x})$ can be
regarded as some smooth functions of position $\mathbf{x}$, the local
linearized Hamiltonian $v_{\mu\nu}\sigma^{\mu}(k^{\nu}-K^{\nu}(\mathbf{x}))$
on the RHS of Eq.(21) can be cast into

$\displaystyle v_{\mu\nu}(\mathbf{x})\sigma^{\mu}(p^{\nu}-\delta
K^{\nu}(\mathbf{x})),$ (22)

where
$\mathbf{p}-\delta\mathbf{K}(\mathbf{x})\equiv\mathbf{k}-\mathbf{K}(\mathbf{x})$,
$\delta\mathbf{K}(\mathbf{x})\equiv\mathbf{K}(\mathbf{x})-\mathbf{K}_{f}$ with
$\mathbf{K}_{f}$ the corresponding Dirac point in strain-free graphene. Unlike
the usual explanation of the strain induced gauge field in grapheneKane ; Neto
, where the effective vector potential is an auxiliary quantity and describes
the mixed effects of both anisotropy of $v_{\mu\nu}$ and the shift of Dirac
point, here the vector potential only represents the relative translation of
the Dirac points, $(e/c)\mathbf{A}(\mathbf{x})=\delta\mathbf{K}(\mathbf{x})$,
its pseudo-magnetic field
$\mathbf{B}(\mathbf{x})=(c/e)\nabla\times\mathbf{K}(\mathbf{x})$, and the
physical effects are mainly determined by the pseudo-magnetic flux through a
loop $(c/e)\oint_{L}\mathbf{K}(\mathbf{x})\cdot d\mathbf{x}$. In the following
sections we shall discuss the properties of $\mathbf{K}(\mathbf{x})$, and
illustrate how a curl field $\mathbf{K}(\mathbf{x})$ is induced by a strain.

From Eq.(16) we know that the vectors representing $t_{1}$,
$\tilde{t}_{2}(\mathbf{K})$, $\tilde{t}_{3}(\mathbf{K})$ in the complex plane
can form a directed triangle for a Dirac point $\mathbf{K}$, as illustrated in
Fig.2a.

Figure 2: (color online). (a) Zero points of $h_{ba}(\mathbf{k})$ determine
two directed triangles with edges $t_{1,2,3}$ in the complex plane. For a
given $t_{1}$ and a fixed direction of $\tilde{t}_{2}$, the arguments of
$\tilde{t}_{3}$ must satisfy conditions (29) to ensure $t_{2,3}\geq 0$. (b)
Dirac points exist if $t_{i}$ satisfy inequalities (23), which describe a
pyramidal domain in $(t_{1},t_{2},t_{3})$ space, if $(t_{1},t_{2},t_{3})$ goes
beyond this domain, an energy gap will be opened.

According to the triangle inequality, we have the following necessary and
sufficient conditions for the existence of the Dirac pointsHasegawa :

$t_{1}+t_{2}\geq t_{3},~{}t_{2}+t_{3}\geq t_{1},~{}t_{3}+t_{1}\geq t_{2}.$
(23)

These conditions define a pyramidal domain in the $(t_{1},t_{2},t_{3})$ space,
shown in Fig.2b. If $t_{1},t_{2},t_{3}$ satisfy inequalities (23), then there
exists two directed triangles with the same edges $t_{1},t_{2},t_{3}$ but
different possible orientations, which determine two angles
$\theta_{1},\theta_{2}$ satisfying
$t_{1}+t_{2}e^{i\theta_{1}}+t_{3}e^{i\theta_{2}}=0$, where
$\theta_{1},\theta_{2}$ are given by the law of cosine

$\displaystyle\begin{split}\theta_{\pm
1}&=\pm\big{[}\pi-\arccos\big{(}\frac{t_{1}^{2}+t_{2}^{2}-t_{3}^{2}}{2t_{1}t_{2}}\big{)}\big{]},\\\
\theta_{\pm
2}&=\pm\big{[}\arccos\big{(}\frac{t_{1}^{2}+t_{3}^{2}-t_{2}^{2}}{2t_{1}t_{3}}\big{)}-\pi\big{]}.\end{split}$
(24)

Thus the Dirac points $\mathbf{K}$ can be determined by letting

$\exp(i\mathbf{K}\cdot\mathbf{a}_{1})=\exp(i\theta_{1}),~{}\exp(i\mathbf{K}\cdot\mathbf{a}_{2})=\exp(i\theta_{2}),$
(25)

so we have

$\mathbf{K}=\frac{1}{2\pi}\big{(}\theta_{1}\mathbf{b}_{1}+\theta_{2}\mathbf{b}_{2}\big{)}+\mathbf{K}_{0},$
(26)

with $\mathbf{b}_{1},\mathbf{b}_{2}$ the reciprocal lattice vectors defined by
$\mathbf{a}_{i}\cdot\mathbf{b}_{j}=2\pi\delta_{ij}$, and
$\mathbf{K}_{0}=n\mathbf{b}_{1}+m\mathbf{b}_{2}$ with $n,m$ are arbitrary
integers. Notice that if $t_{1}+\tilde{t}_{2}+\tilde{t}_{3}=0$, then
$t_{1}+\tilde{t}^{*}_{2}+\tilde{t}^{*}_{3}=0$, this implies that there exists
two Dirac points $\mathbf{K}(\mathbf{x})$ and $-\mathbf{K}(\mathbf{x})$.
However, if $(t_{1},t_{2},t_{3})$ exactly locates on the boundary surface of
the pyramid, e.g., $t_{1}=t_{2}+t_{3}$, then the two triangles will mutually
coincide and $\tilde{t}_{2}=\tilde{t}^{*}_{2}$,
$\tilde{t}_{3}=\tilde{t}^{*}_{3}$(see Fig.2a), hence $\mathbf{K}(\mathbf{x})$
and $-\mathbf{K}(\mathbf{x})$ become equivalent, and
$\vec{\alpha}=i(t_{2}\mathbf{a}_{1}+t_{3}\mathbf{a}_{2})$ becomes a pure
imaginary vector, so the Fermi velocity in the directions perpendicular to
$\vec{\alpha}$ vanishes(Fig.3b and 3d)Montambaux08 ; Montambaux09 ; Pereira .
If $(t_{1},t_{2},t_{3})$ goes beyond the domain defined by Eq.(23), e.g.,
$t_{1}>t_{2}+t_{3}$, Eq.(16) will have no any root, an energy gap with
magnitude $E_{g}=2(t_{1}-t_{2}-t_{3})$ will occur at the corresponding points
$\mathbf{K}_{\pm}=\pm 1/2(\mathbf{b}_{1}+\mathbf{b}_{2})$(see Fig.3c)Zhou ,
and the effective Hamiltonian (22) must be further modified by adding a mass
term and some second order terms.

Figure 3: (color online) (a) Energy band when two Dirac points are very close,
where $t_{1}=2.8,t_{2,3}=1.45,$ (b), (d) $t_{1}=2.8,t_{2,3}=1.4,$ two Dirac
points are equivalent(superposed), (c) $t_{1}=2.8,t_{2,3}=1.35,$ an energy gap
occurs.

Another important property is the range of $\mathbf{K}(\mathbf{x})$. We shall
prove that the Dirac points must be confined within some special regions of
the BZ. To this end, notice that if a Dirac point
$\mathbf{K}=(1/2\pi)(\theta_{1},\theta_{2})$ is given, then its associated
$t_{1,2,3}$ can also be determined up to an arbitrary factor, except six
special cases of $\theta_{1},\theta_{2}=0,\pm\pi$(see Fig. 2a). If
$\theta_{1},\theta_{2}\neq 0,\pm\pi$. According to the law of sines we have

$\displaystyle\frac{t_{2}}{t_{1}}=\frac{\sin\theta_{2}}{\sin(\theta_{1}-\theta_{2})},~{}~{}\frac{t_{3}}{t_{1}}=\frac{\sin\theta_{1}}{\sin(\theta_{2}-\theta_{1})},$
(27)

or

$\displaystyle(t_{1},t_{2},t_{3})\propto(\sin(\theta_{2}-\theta_{1}),-\sin\theta_{2},\sin\theta_{1}).$
(28)

For the six special cases we have: if $(\theta_{1},\theta_{2})=\pm(\pi,\pi),$
$t_{1}=t_{2}+t_{3};$ if $(\theta_{1},\theta_{2})=\pm(\pi,0),$
$t_{2}=t_{1}+t_{3};$ if $(\theta_{1},\theta_{2})=\pm(0,\pi),$
$t_{3}=t_{1}+t_{2}$. From Eq.(28) we can find that the $\theta_{1},\theta_{2}$
must satisfy some constrain conditions to guarantee $t_{1,2,3}\geq 0$, as
illustrated in Fig.2a. For an arbitrary $t_{2}$ and a fixed
$\theta_{1}$(direction of $\tilde{t}_{2}$), $\tilde{t}_{3}$ must point in a
direction between the directions of $-\tilde{t}_{2}$ and negative real axis,
i.e., argument $\theta_{1},\theta_{2}$ must satisfy

$\displaystyle\begin{split}\theta_{1}&+\pi<\theta_{2}<\pi,&\theta_{1}&\in(-\pi,0),\\\
-\pi&<\theta_{2}<\theta_{1}-\pi,&\theta_{1}&\in(0,\pi).\\\ \end{split}$ (29)

These two inequalities respectively determine the range of
$\mathbf{K}(\mathbf{x})$ and $-\mathbf{K}(\mathbf{x})$. They describe two open
triangles $\bigtriangleup MM^{\prime}_{1}M^{\prime\prime}_{1}$ and
$\bigtriangleup M^{\prime}M_{1}M^{\prime\prime}_{3}$ in reciprocal space, as
shown in Fig.4, each one has 1/8 the area of a unit cell of the reciprocal
space(the parallelogram
$M^{\prime\prime}M^{\prime\prime}_{1}M^{\prime\prime}_{2}M^{\prime\prime}_{3}$),
and each Dirac point is confined within a triangle, so, the Dirac points
$\mathbf{K}$ and $-\mathbf{K}$ can meet(become equivalent) only at the
vertexes of $\bigtriangleup MM^{\prime}_{1}M^{\prime\prime}_{1}$ and
$\bigtriangleup M^{\prime}M_{1}M^{\prime\prime}_{3}$. The remaining
hexagon(blue in Fig.4) is a forbidden region for the Dirac points. This
confinement also limits the order of magnitude of
$\nabla\times\mathbf{K}(\mathbf{x})$, i.e., the strain induced pseudo-magnetic
field.

Figure 4: (color inline) Rang of the Dirac points consists of six triangles in
first the BZ, or $\bigtriangleup MM^{\prime}_{1}M^{\prime\prime}_{1}$ and
$\bigtriangleup M^{\prime}M_{1}M^{\prime\prime}_{3}$.

In order to show the underlying regularity, here we have ignored the
variations of $\mathbf{b}_{1},\mathbf{b}_{2}$ with the deformation of lattice,
and simply sketch all $\mathbf{K}=(k_{1},k_{2})$ in the same affine frame.
After translating to the first BZ of graphene, $\bigtriangleup
MM^{\prime}_{1}M^{\prime\prime}_{1}$ and $\bigtriangleup
M^{\prime}M_{1}M^{\prime\prime}_{3}$ are equivalent to a ringlike region
consists of six triangles $\triangle MKM^{\prime}_{1},$ $\triangle
M_{1}K_{2}M^{\prime\prime}_{2},$ $\triangle M^{\prime\prime}K_{1}M^{\prime},$
etc.

In order to illustrate how a non-vanishing $\nabla\times\mathbf{K}$ is induced
by strain, we only need to analyze three ideal cases, in which only one
$t_{i}$ is slightly changed, $t_{i}\rightarrow t_{0}+\delta t_{i}$, while the
other two $t_{j,k}$ remain constant, $t_{j}=t_{k}=t_{0}$, which can also be
roughly regarded as that the bond $\mathbf{c}_{i}$ is elongated(or compressed)
while the other two bonds $\mathbf{c}_{j},\mathbf{c}_{k}$ and their directions
remain fixed(see Fig.5b). Notice that the Dirac points only depend on the
relative proportions of $t_{1},t_{2},t_{3}$, so, as an equivalent case, we can
always assume that $t_{1}$ remains constant and only $t_{2}$, $t_{3}$ are
variables. Moreover, in these equivalent cases the $\tilde{t}_{2}$ and
$\tilde{t}_{3}$ can be determined by the end of the vector
$t_{1}+\tilde{t}_{2}=-\tilde{t}_{3}$, denoted by $P$ in Fig.5a. So, we can
represent the variation of the Dirac points by the shift of the point P. To
this end, we have to determine the corresponding P of the three classes of
characteristic points in the range of the Dirac points: (1) $K$(or
$K^{\prime}$) etc.(see Fig.4), according to Eq.(27), Dirac points locate at
these two points only if $t_{1}=t_{2}=t_{3}$, their corresponding P is located
at $K$(or $K^{\prime}$) in Fig.5a; (2) critical points $M(M_{1})=(\pm 1/2,0)$,
$M^{\prime}(M^{\prime}_{1})=(0,\pm 1/2)$, and
$M^{\prime\prime}(M^{\prime\prime}_{2})=(\pm 1/2,\pm 1/2),$ in these cases
there exist only one Dirac point since the points $\mathbf{K}$ and
$-\mathbf{K}$ are equivalent, their corresponding $(t_{1},t_{2},t_{3})$ are
located on the boundary of the pyramidal domain, while their corresponding P
are located at the real axis in Fig.5a; (3) $O$, $O^{\prime}$,
$O^{\prime\prime}$ etc., their corresponding P are the centers of two circles
and the infinite limit points of the straight line $KK^{\prime}$ in Fig.5a,
which respectively correspond to the limits of $t_{2}\rightarrow 0$,
$t_{3}\rightarrow 0$ and $t_{1}\rightarrow 0$(equivalent to
$t_{2}=t_{3}\rightarrow\infty$).

Figure 5: (a) The trajectories of P in three ideal cases. (b) A schematic
diagram of the shift of Dirac point $\delta\mathbf{K}$ and its curl,
$\delta\mathbf{K}$ perpendicular to $\mathbf{c}_{i}$, if bond $\mathbf{c}_{i}$
is slightly elongated.

Now we analyze the shifts of the Dirac points in the three ideal situations.
(1) $t_{1},t_{2}$ remain constant and $t_{2}=t_{1}$, only $t_{3}$ is variable,
the trajectory of the corresponding P is a circle with radius $t_{1}$ and
centered at the point $(t_{1},0)$ in Fig.5a, so the arguments
$(\theta_{1},\theta_{2})$ satisfy

$\displaystyle\begin{split}\theta_{1}-2\theta_{2}-2\pi=0,&~{}~{}\theta_{1}\in(-\pi,0),\\\
\theta_{1}-2\theta_{2}+2\pi=0,&~{}~{}\theta_{1}\in(0,\pi).\end{split}$ (30)

They describe line segments $M^{\prime}O_{1}$ and $M^{\prime}_{1}O$ in Fig.4;
(2) $t_{3}(=t_{1})$ remain constant while $t_{2}$ is variable, the trajectory
of corresponding P is another circle with radius $t_{1}$ centered at the
origin, its associated $(\theta_{1},\theta_{2})$ satisfy

$\displaystyle\begin{split}\theta_{2}-2\theta_{1}+2\pi=0,&~{}~{}\theta_{2}\in(-\pi,0),\\\
\theta_{2}-2\theta_{1}-2\pi=0,&~{}~{}\theta_{2}\in(0,\pi),\end{split}$ (31)

which describe $MO^{\prime}$ and $M_{1}O^{\prime}_{1}$ in Fig.4; (3) $t_{1}$
remains constant while $t_{2},t_{3}$ are variable but $t_{2}=t_{3}$(or vice
versa, $t_{1}$ is variable, $t_{2}(=t_{3})$ remain constant), the trajectory
of P is straight line $KK^{\prime}$, $(\theta_{1},\theta_{2})$ satisfy

$\theta_{1}+\theta_{2}=0,~{}~{}\frac{\pi}{2}<|\theta_{1}|<\pi,$ (32)

which describe $M^{\prime\prime}_{1}O^{\prime\prime}$ and
$M^{\prime\prime}_{3}O^{\prime\prime}_{1}$ in Fig.4. Summarizing
Eqs.(30)(31)(32) and comparing with Fig.4, we observe that if a band, e.g.,
$\mathbf{c}_{1}$ is slightly elongated (or compressed) along its direction,
$\mathbf{c}_{1}\rightarrow(1+\delta)\mathbf{c}_{1}$, while the other two bonds
$\mathbf{c}_{2},\mathbf{c}_{3}$ remain fixed, then $t_{1}$ will be slightly
changed while $t_{2}(=t_{3})$ remain constant, the Dirac point $\mathbf{K}$
will be slightly moved in the direction perpendicular to $\mathbf{c}_{1}$,
i.e., $\delta K_{y}\neq 0$(see Fig.4, $K$ moves towards $O^{\prime\prime}$ if
$t_{1}$ decreases, towards $M^{\prime\prime}_{1}$ if $t_{1}$ increases),
$\mathbf{K^{\prime}}(=-\mathbf{K})$ is moved in the opposite direction. Thus,
if the elongation of $\mathbf{c}_{1}$ is slowly varying in the x-direction,
i.e, $\partial t_{1}/\partial x\neq 0$, then $\partial K_{y}/\partial x\neq
0$, the other two cases are similar. So, a nonuniform strain as schematically
shown in Fig.5b can induce a curl field $\mathbf{K}(\mathbf{x})$,
$\nabla\times\mathbf{K}\neq 0$.

###### Acknowledgements.

We thank Yugui Yao, Chengshi Liu, Yichun Liu for their helpful discussions.
This work was supported by the National Science Foundation of China(Grant
Nos.10974027, 50725205, 50832001).

## References

  * (1) K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, Y. Zhang, S. V. Dubonos, I. V. Gregorieva, and A. A. Firsov, Science 306, 666(2004); K. S. Novoselov, D. Jiang, F. Schedin, T. J. Booth, V. V. Khotkevich, S. M. Morozov, A. K. Geim, Proc. Natl. Acad. Sci. 102, 10451(2005).
  * (2) M. A. H. Vozmediano, M. P. Lopez-Sancho, and F. Guinea, Phys. Rev. Lett. 89, 166401(2002).
  * (3) J.C. Meyer, K. Geim, M.I. Katsnelson, K.S. Novoselov, T.J. Booth, and S. Roth, Nature 446, 60(2007).
  * (4) A.K. Geim, and K.S. Novoselov, Nature Materials 6, 183(2007); C.W.J. Beenakker, Rev. Mod. Phys. 80, 1337(2008).
  * (5) A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, and A. K. Geim, Rev. Mod. Phys. 81, 109(2009).
  * (6) G. W. Semenoff, Phys. Rev. Lett. 53, 2449(1984).
  * (7) V. P. Gusynin, and S. G. Sharapov, Phys. Rev. Lett. 95, 146801(2005); V. P. Gusynin, V. A. Miransky, and S. G. Sharapov, Phys. Rev. B 74, 195429(2006).
  * (8) K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, M. I. Katsnelson, I. V. Gregorieva, S. V. Dubonos, and A. A. Firsov, Nature 438, 197(2005).
  * (9) Y. Zhang, Y.-W. Tan, H. L. Stormer, and P. Kim,Nature 438, 201(2005).
  * (10) F. Guinea, M. I. Katsnelson, and A. K. Geim, Nuatur Physics 6, 30(2010).
  * (11) N. Levy, S.A.Burke, K.L.Meaker, M.Panlassigui, A.Zettl, F.Guinea, and A.H. Crommie, Science 329, 544 (2010).
  * (12) S. L. Zhu, B. Wang, and L.-M. Duan, Phys. Rev. Lett. 98, 260402(2007); O. Bahat-Treidel, O. Peleg, M. Grobman, N. Shapira, M. Segev, and T. Pereg-Barnea, Phys. Rev. Lett. 104, 063901(2010).
  * (13) Y. Hasegawa, R. Konno, H. Nakano, and M. Kohmoto, Phys. Rev B 74, 033413(2006).
  * (14) P. Dietl, F. Piéchon, and G. Montambaux, Phys. Rev. Lett. 100, 236405(2008).
  * (15) G. Montambaux, F. Pi chon, J.-N. Fuchs, and M. O. Goerbig, Phys. Rec. B 80, 153412(2009); G. Montambaux, F. Pi chon, J.-N. Fuchs and M. O. Goerbig, Eur. Phys. J. B 72, 509(2009).
  * (16) P. R. Wallace, Phys. Rev. 71, 622(1947).
  * (17) C.L.Kane and E.J.Mele Phys. Rev. Lett. 78, 1932(1997); H. Suzuura and T. Ando, Phys. Rev. B 65, 235412(2002); J.L.Mañe, Phys. Rev. B 76, 045430(2007).
  * (18) V.M. Pereira, A.H. Castro Neto, and N.M.R. Peres, Phys. Rev. B 80, 045401(2009).
  * (19) Other mechanisms of gap opening see e.g., S.Y. Zhou, G.H.Gweon, A.V. Fedorov, P.N.First, W.A.De Heer, D.H.Lee, F.Guinea, A.H.Castro Neto, and A.Lanzara, Nature Materials 6, 770(2007); G. Giovannetti, P. A. Khomyakov, G. Brocks, P.J. Kelly, and J. van den Brink, Phys. Rev. B 76, 073103(2007); R. Martinazzo, S. Casolo, and G. F. Tantardini, Phys. Rev. B 81, 245420(2010).


figurec

# A preliminary assessment of the sensitivity of uniaxially-driven fusion
targets to flux-limited thermal conduction modeling

D. A. Chapman<EMAIL_ADDRESS>First Light Fusion Ltd., Unit
9/10 Oxford Industrial Park, Mead Road, Yarnton, Kidlington OX5 1QU, United
Kingdom  J. D. Pecover  First Light Fusion Ltd., Unit 9/10 Oxford Industrial
Park, Mead Road, Yarnton, Kidlington OX5 1QU, United Kingdom  N. Chaturvedi
First Light Fusion Ltd., Unit 9/10 Oxford Industrial Park, Mead Road, Yarnton,
Kidlington OX5 1QU, United Kingdom  Centre for Inertial Fusion Studies,
Blackett Laboratory, Imperial College, London SW7 2AZ, United Kingdom  N.
Niasse  First Light Fusion Ltd., Unit 9/10 Oxford Industrial Park, Mead Road,
Yarnton, Kidlington OX5 1QU, United Kingdom  M. P. Read  First Light Fusion
Ltd., Unit 9/10 Oxford Industrial Park, Mead Road, Yarnton, Kidlington OX5
1QU, United Kingdom  D. H. Vassilev  First Light Fusion Ltd., Unit 9/10 Oxford
Industrial Park, Mead Road, Yarnton, Kidlington OX5 1QU, United Kingdom  J. P.
Chittenden  Centre for Inertial Fusion Studies, Blackett Laboratory, Imperial
College, London SW7 2AZ, United Kingdom  N. Hawker  First Light Fusion Ltd.,
Unit 9/10 Oxford Industrial Park, Mead Road, Yarnton, Kidlington OX5 1QU,
United Kingdom  N. Joiner  First Light Fusion Ltd., Unit 9/10 Oxford
Industrial Park, Mead Road, Yarnton, Kidlington OX5 1QU, United Kingdom

(August 28, 2024)

###### Abstract

The role of flux-limited thermal conduction on the fusion performance of the
uniaxially-driven targets studied by Derentowicz et al.
Derentowicz_JourTechPhys_1977  is explored as part of a wider effort to
understand and quantify uncertainties in ICF systems sharing similarities with
First Light Fusion’s projectile-driven concept. We examine the role of
uncertainties in plasma microphysics and different choices for the numerical
implementation of the conduction operator on simple metrics encapsulating the
target performance. The results indicate that choices which affect the
description of ionic heat flow between the heated fusion fuel and the gold
anvil used to contain it are the most important. The electronic contribution
is found to be robustly described by local diffusion. The sensitivities found
suggest a prevalent role for quasi-nonlocal ionic transport, especially in the
treatment of conduction across material interfaces with strong gradients in
temperature and conductivity. We note that none of the simulations produce
neutron yields which substantiate those reported in
Ref.Derentowicz_JourTechPhys_1977 , leaving open future studies aimed at more
fully understanding this class of ICF systems.

##  I Introduction

First Light Fusion (FLF) is investigating a novel approach to controlled
inertial confinement fusion (ICF) using efficient, low-cost, projectile-based
driver technology. A key characteristic of the concept is that the implosion
process initially proceeds with a single, strong shock directed along a single
axis, although the dynamics later-in-time can have convergent, non-planar
aspects. As with all routes to ICF, a crucial element of target design and
optimization is the validation of the numerical tools, and the multitude of
options which underpin them. Sensitivity studies seek to assess the impact of
factors such as code configuration choices, numerical methods and the basic
properties of materials and are crucial for enhancing the confidence in
simulation-led predictions. Such studies are becoming commonplace in the ICF
community Gaffney_NuclFusion_2013 ; Gaffney_HEDP_2013 ;
Melvin_PhysPlasmas_2015 ; Guymer_PhysPlasmas_2015 .

For the planar driver geometry of interest there are very few experiments
accessing the fusion regime against which code benchmarking can be undertaken.
One notable class of experiments where the fuel is collapsed directly by a
single planar shock are are those of Derentowicz et al.
Derentowicz_JourTechPhys_1977 , which feature a fuel-filled conical cavity in
a metal anvil driven by a uniform, planar Mach wave produced by the implosion
of a conical liner. Several variations of this experiment have been performed
Krasyuk_QuantElectronics_2005  that consider different driver technologies,
such as electrical discharge explosion Shyam_AtomkernenergieKerntechnik_1984 ,
direct laser ablation Vovchenko_JETPLett_1977 ; Mason_ApplPhysLett_1979  and
relativistic electron beams Bogolyubskii_JETPLett_1976 . The reported output
of these targets of $10^{4}-10^{7}$ neutrons is of substantial utility from
the perspective of validating both integrated simulations and diagnostic
calibration. Of these, the original concept provides the simplest and
therefore most valuable validation case.

As for an ICF experiment, the modeling is challenging: compressible
hydrodynamics, thermal conduction, viscous drag, radiation generation,
transport and absorption and detailed thermophysical material properties
across a wide range of density-temperature space must all be considered.
Initial examinations have shown that the most important process for
determining fuel energetics is thermal conduction, which acts to homogenize
strongly-localised heating resulting from reverberating shock waves.

For plasmas featuring gentle thermal gradients the rate of conduction loss
from a small volume of heated fuel is dictated principally by the thermal
conductivities of the electrons and ions; $\kappa_{e}$ and $\kappa_{i}$,
respectively. In low-density, high-temperature, hydrogen-like systems the
electronic contribution dominates due to their higher mobility, with ionic
conduction only competing when $T_{i}\gtrsim 4T_{e}$ Larsen_book . Since
shock-heating couples most efficiently to the ions, nonequilibrium conditions
where both electronic and ionic conduction are important are readily
attainable. Additionally, electron-ion temperature relaxation directly affects
not only the plasma reactivity but also influences the balance between
electronic and ionic conduction loss. Although the theoretical descriptions of
these properties are well-known under the near-ideal conditions of the fuel,
the solid target components are predicted to be driven into warm dense matter
(WDM) states, with close-to-solid densities and temperatures of a few
$\mathrm{eV}$. Under such extreme conditions even basic thermophysical
properties can be substantially uncertain Graziani_book .

The system studied in Ref. Derentowicz_JourTechPhys_1977  uses a planar shock
in copper to drive the collapse of a deuterium-filled conical cavity with a
full internal angle at the tip of $60{}^{\circ}$. Between the copper and the
fuel is a polythene layer, referred to in this work as the ‘coverslip’. The
material in which the cavity is made is gold and is referred to here as the
‘anvil’. At the stagnation point of the system, where the converging flow into
the tip of the cavity is reflected from the anvil, a transient nonequilibrium
state is produced which sets up a thermal wave across the interface between
the fuel and interior anvil wall. The resulting heat flow is driven by both
electronic and ionic conduction and occurs in the presence of sharp density
and temperature gradients between regions with very different conductivities
and heat capacities. Accurately capturing the heat flux under such conditions
is in general a problem which can only be fully understood using state-of-the-
art kinetic codes (see, e.g. Refs. Larroche_EuroPhysJourD_2003 ;
Larroche_PhysPlasmas_2012 ; Taitano_JourComputPhys_2018 ) and is one of the
many challenges of present concern in the ICF community
Rinderknecht_PlasmaPhysControlFusion_2018 . Reduced kinetic models, such as
the Schurtz-Nicolaï-Busquet (SNB) model Schurtz_PhysPlasmas_2000 ;
Cao_PhysPlasmas_2015  are available which reproduce many of the features
expected from heat flow in the kinetic regime Hoffman_PhysPlasmas_2015 ;
Brodrick_PhysPlasmas_2017 ; Sherlock_PhysPlasmas_2017 , although recent
experiments Henchen_PhysRevLett_2018 ; Henchen_PhysPlasmas_2019  have
demonstrated even larger heat flow suppression than predicted.

Commonly-used simple models, such as local flux-limited thermal conduction
Bell_PhysRevLett_1981 , offer a computationally tractable method used in many
ICF simulation codes Meezan_PhysPlasmas_2020 . Although conceptually
straightforward, the precise manner in which they are implemented is subject
to uncertainty due to the large number of free parameters involved and
adequate evidence-based values are seldom available. Of particular importance
are: the choices of the magnitude of the flux-limited heat flow, i.e. the
values of the flux limiter coefficients; the manner in which the transition
between purely diffusive and flux-limited conduction is undertaken; and how
the relevant quantities are to be utilized within the numerical framework of
the host code, e.g. in the spatial discretization scheme.

In this work we report on the initial results of an ongoing simulations-based
study investigating the influence of uncertainties in FLF’s predictive
modeling capability using the experimental results from Ref.
Derentowicz_JourTechPhys_1977  as a benchmark. The scope of this paper is
focused on the impact of flux-limited thermal conduction on an idealized model
of the target, as this is currently believed to be the largest effect on the
fusion output. With suitably well-converged simulations we establish a
reference case against which the impact of configuration changes can be
quantifiably measured using a few simple metrics. The basic form of the
conduction operator is described, with emphasis given to the main sources of
uncertainty arising from free choices in its implementation. The changes with
the strongest influence on the evolution of the target dynamics and fuel
energetics are found to be the accuracy of the electron thermal conductivity
of the fuel, the magnitude of the maximum ionic heat flux and the treatment of
interface conduction, both of which can only be fully understood with more
advanced tools. Since the thermophysical properties of the fuel under the
states produced are thought to be well-known, we conclude that an improved
understanding of ionic conduction in the kinetic regime is therefore crucial
if we are to be able to use this experimental platform as a means to validate
simulations of relevance to FLF’s mission.

##  II Simulation setup

###  II.1 Idealized target model

(a)

(b)

(c)

Figure 1:  (a): Simulation geometry and EGrid zoning for the idealized model,
featuring a shocked Cu driver, PMMA coverslip, deuterium fuel and gold anvil.
The left-hand part of the plot shows a logarithmic map of mass density whilst
the right-hand side shows the material ‘color’ indicator. The tip of the
conical cavity is rounded with a radius of curvature of $10\,\mathrm{\mu m}$,
as shown in detail in the inset. (b): Time histories of the BWA ion (upper
panel) and electron (middle panel) temperatures and mass density (lower
panel). The regions with the gray background indicate the time window over
which simulation convergence is assessed. (c): Convergence metric
$\epsilon_{\text{tot}}$ as a function of minimum cell size.

The present work is undertaken using Hytrac; one of the two in-house
radiation-hydrodynamics codes developed by FLF. The code is based on the front
tracking approach Jian_JComputPhys_2006  and also implements cell-based
adaptive mesh refinement, similar to, e.g., RAGE Gittings_ComputSciDisc_2008 .
The multi-physics model includes two-temperature conduction, viscous fluxes
Vold_PhysPlasmas_2017  and radiation transport via the $P_{1/3}$-AFL method
Mihalas_JQSRT_1982 ; McGlinchey_thesis_2017 . The code further supports either
analytic or tabulated models for material equation of state (EoS), plasma
microphysics and fusion reactivity.

The work presented in Ref. Derentowicz_JourTechPhys_1977  shows that the drive
into the fuel is well-approximated with a single, uniform shock into the
plastic coverslip. The pressure in the copper is inferred from shock velocity
measurements to be $46\,\mathrm{Mbar}$, which leads to a shock in coverslip
with a pressure of $13\,\mathrm{Mbar}$. The latter is modeled as PMMA as a
surrogate for polythene due to difficulties matching the principal Hugoniot
using the FEOS EoS model Faik_ComputPhysCommun_2018 . Since the densities of
the two materials are sufficiently similar this compromise is believed to be
reasonable.

The target mesh geometry used for our simulations is shown in Fig. 1a. The
number of AMR refinement levels is 4, giving a minimum cell size of
$1.7\,\mathrm{\mu m}$ on the Eulerian grid (EGrid). The inset shows how the
cavity tip has had a radius of curvature of $10\,\mathrm{\mu m}$ added to
account for expected manufacturing tolerances. The shocked copper region,
denoted as the ‘driver’, at the top of the domain is initialised from the
thermodynamic state resulting from shocking initially ambient copper to the
experimentally inferred pressure along the principal Hugoniot. All other
materials are initialised at a pressure of $1\,\mathrm{bar}$ and their
appropriate ambient mass densities to produce initially static interfaces.

In these simulations the viscous flux and radiation transport operators are
not used. In the case of viscous effects this is because the predominant
influence is on the damping of the incident shock and not the energetics of
the cavity collapse, from which the bulk of the yield is expected. Full
radiation-hydrodynamics modelling suggests the internal cavity dynamics may be
influenced by bremsstrahlung emission from the fuel to a small degree through
radiative ablation of the wall, although the conduction losses are vastly more
important to the energetics. Disabling these operators also makes tractable a
far larger study focused on the physics of conduction, which is our principal
focus here. More details of the cavity collapse dynamics and the impact of the
full range of multi-physics options will be discussed in future work.

###  II.2 Assessment of convergence

Assessment of the impact of changes to the setup of the multi-physics model
must be done on the basis that any differences observed can confidently be
attributed to the predicted influence of the options, rather than due to lack
of simulation convergence. We pursue convergence through variation of the
EGrid base resolution at a constant refinement level and comparing the
difference between key metrics of progressively higher-resolution runs.

The metrics used for this purpose are the burn-weighted average (BWA)
histories of the ion and electron temperatures and mass density

$\displaystyle\\!\left\langle X\right\rangle_{b}\\!\left(t\right)=$
$\displaystyle\,\dfrac{\displaystyle\int_{V}\mathrm{d}\mathbf{r}\,\,\dfrac{\partial^{2}Y_{\text{n}}}{\partial
V\partial
t}\,X(\mathbf{r},t)}{\displaystyle\int_{V}\mathrm{d}\mathbf{r}\,\,\dfrac{\partial^{2}Y_{\text{n}}}{\partial
V\partial t}}\,.$ (1)

Here, $X$ stands for the quantity of interest,
$\partial^{2}Y_{\text{n}}/\partial V\partial t$ represents the total neutron
yield emission per unit volume per unit time and $V$ denotes the volume
occupied by the fusion fuel. The convergence parameter $\epsilon$ with respect
to variable $X$ for a simulation with minimum cell size $\Delta_{\min}$
relative to a suitable reference case is then defined as

$\displaystyle`\epsilon_{X}=$
$\displaystyle\,\\!\left[\frac{\sum_{i_{\min}}^{i_{\max}}\\!\left(\\!\left\langle
X\right\rangle_{b}\\!\left(t_{i};\Delta_{\min}\right)-\\!\left\langle
X\right\rangle_{b}\\!\left(t_{i};\Delta_{\min}^{\text{ref}}\right)\right)^{2}}{i_{\max}-i_{\min}}\right]^{1/2}\,.$
(2)

Here, the index $i$ runs between times bounded by the onset ($t_{i_{\min}}$)
and secession ($t_{i_{\max}}$) of neutron production. These are respectively
defined as the earliest and latest times where the neutron production rate
exceeds $10\%$ of its peak value. The convergence metric is then defined as
the quadratic average over all variables, i.e.
$\epsilon_{\text{tot}}=\\!\left(\sum_{X}\epsilon_{X}^{2}/\sum_{X}\right)^{1/2}$.
For pure deuterium fuel at the modest temperatures predicted one may
approximately write

$\displaystyle\frac{\partial^{2}Y_{\text{n}}}{\partial V\partial t}\approx$
$\displaystyle\,\frac{n_{\text{D}}^{2}}{2}\\!\left\langle\sigma
v\right\rangle_{d(D,{}^{3}\text{He})\text{n}}\,.$ (3)

In Eq. (3) the thermal reactivity $\\!\left\langle\sigma v\right\rangle$ is
given by fits due to Bosch and Hale Bosch_NuclFusion_1992 .

A strong convergence trend can be seen in both the timing and shape of the
peaks in $\\!\left\langle T_{i}\right\rangle_{b}$ and $\\!\left\langle
T_{e}\right\rangle_{b}$ with increasing base resolution (Fig. 1b), indicating
that the hydrodynamics of the cavity implosion are properly captured in our
simulations. This is further exemplified in Fig. 1c. The two main spikes in
the BWA ion temperature correspond to the reflection of the incident shock
from the wall and tip of the cavity. The subsequent reverberations in the
collapsing cavity produce a long-lived plasma in which conduction losses are
roughly balanced by the heating due to compression. The electrons do not show
the same strongly peaked structure since they are treated as incompressible
and are therefore isentropically heated through the shock reflection.
Nonequilibrium states with $T_{i}/T_{e}>3$ are achieved in this period before
equilibration to a temperature of $T\sim 160-180\,\mathrm{eV}$, which persists
for around $10\,\mathrm{ns}$. The role of ion conduction in dissipating the
energy imparted to the fuel at the tip of the cavity is therefore likely to be
important, as is the rate of temperature equilibration. Overall, the strength
of the observed trends indicates that our results at a base resolution of 40
cells (smallest cell size of $1.7\,\mathrm{\mu m}$) provides a suitable
reference case for quantifying the impact of changes due to conduction
physics, without conflation with issues related to lack of convergence.

##  III Conduction modeling in ICF simulations

###  III.1 Flux-limited thermal conduction

It is well-known that the modeling of losses due to thermal conduction
constitute a crucial aspect of the power balance in prospective ICF systems
Atzeni_book . A commonly-used implementation strategy for the conduction
operator in radiation-hydrodynamics codes follows the local flux-limited
diffusion approach to Fourier’s law

$\displaystyle\mathbf{Q}_{\text{Four.}}=$ $\displaystyle\,-\kappa\nabla T\,,$
(4)

the magnitude of which is restricted to a fraction, $\alpha$, of the free-
streaming limit

$\displaystyle Q_{\max}=\alpha
n\,v_{\text{th}}\,k_{\text{B}}T;\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ v_{\text{th}}=(k_{\text{B}}T/m)^{1/2}\,,$ (5)

in order to prevent unphysically large fluxes from arising in regions with
steep temperature gradients. Such conditions arise routinely near the ablation
front of laser-irradiated solids, e.g. in short-pulse laser-driven targets
Brown_PhysRevLett_2011 ; Hoarty_PhysRevLett_2013  or the heating of hohlraum
walls at NIF Meezan_PhysPlasmas_2020 . In the present study sharp temperature
gradients are expected to occur in proximity to strong shocks and also at
interfaces between materials with very different heat capacities.

The flux-limited model represents a phenomenological correction to purely
diffusive (local) energy transport and does not capture the preheating effect
arising from long-range (near collisionless) propagation of the high-energy
tail of the particle distribution. The impact of nonlocal transport is now of
growing concern in the ICF community Hoffman_PhysPlasmas_2015 ;
Rinderknecht_PlasmaPhysControlFusion_2018  and has driven the development of
numerous high-fidelity kinetic simulation tools Larroche_EuroPhysJourD_2003 ;
Huang_PhysPlasmas_2017  and accurate reduced models Schurtz_PhysPlasmas_2000 .
Unfortunately, these more accurate capabilities are seldom practical for
implementation in fully-integrated simulations and the flux-limited approach
remains in common usage. Even for this simplest model, however, there are a
number of free parameters to set and numerical choices to be made;
independently for the electrons and ions.

###  III.2 Flux limiter coefficients

For the electrons, the value of $\alpha_{e}$ is typically tuned to match
experimental data or designed to match results from high-fidelity simulations,
such as kinetic codes, with $\alpha_{e}\approx 0.05$ being a common choice
informed by early Vlasov-Fokker-Planck (VFP) results Bell_PhysRevLett_1981 .
Despite being informed by numerical results a broad range of values may be
needed to match different experimental data sets, with values up to
$\alpha_{e}=0.15$ (the so-called ‘high-flux’ model Rosen_HEDP_2011 ) being
required. On the other hand, recent work by Meezan et al.
Meezan_PhysPlasmas_2020  has shown that a flux-limited approach with a lower
value of $\alpha_{e}=0.03$ can qualitatively explain results from several
independent diagnostics simultaneously. Thus, a wide range of values for this
parameter can be substantiated for sensitivity studies.

The ion heat flow is often either not flux-limited or its contribution is
neglected entirely on the basis that electronic heat flow is usually dominant.
The value of $\alpha_{i}$ is generally taken to be much larger than
$\alpha_{e}$, although there is presently little-to-no consensus on a
canonical value. Indeed, ionic VFP simulations strongly discredit the idea
that a single value of $\alpha_{i}$ can universally describe ionic heat flow,
especially in converging shocks Larroche_EuroPhysJourD_2003 . The default
value used in this work is taken to be $0.5$ as this has been found in several
other radiation-hydrodynamics codes, such as HYADES Larsen_JQSRT_1994  and
DEIRA DEIRA_code_manual . We consider the value of $\alpha_{i}$ to be one of
the most uncertain aspects of our modeling.

###  III.3 Interpolation functions

Another important free choice in the implementation of flux-limited conduction
is the manner in which the capping of the Fourier heat flow to the maximum
value given by Eq. (5) is achieved. A particularly prevalent choice in the
literature is to simply take the minimum value Atzeni_book ; DEIRA_code_manual

$\displaystyle Q_{\text{eff}}=$
$\displaystyle\,\min\\!\left(Q_{\text{Four.}},Q_{\text{max}}\right)\,,$ (6)

which leaves the Fourier heat flow unmodified until the transition point. The
impact of such a hard switch is discussed by Meezan et
al.Meezan_PhysPlasmas_2020 , wherein the propagation of a nonlinear thermal
wave is studied for $\alpha_{e}=0.03-0.25$. The sharp front found for smaller
flux-limiter coefficients results directly from the discontinuous profile in
the heat flux due to Eq. (6), which is unlikely to be representative of the
true solution. Indeed, our simple conduction verification problems have shown
that this approach can lead to spurious break-up of the wave front and often
shows undesirable numerical noise (see the related curve in Fig. 2a).

(a)

(b)

(c)

Figure 2:  (a): Spatial temperature profile for a simple nonlinear 1D thermal
wave propagation test in an ideal gas with Spitzer-like thermal
conductivitynrl_formulary , comparing different options for the interpolation
functions (b) and for evaluating quantities on the cell faces (c). The
harmonic mean case (dashed pink curve) moves just over $2.5\,\mathrm{\mu m}$
from the initial temperature discontinuity ($x_{0}=5\,\mathrm{\mu m}$) and
maintains a temperature of over $250\,\mathrm{eV}$ behind the front. (b):
Interpolation functions $f_{\text{interp}}$ governing the transition between
$\kappa_{\text{tab}}$ and $\kappa_{\max}$. (c): Interpolation functions for
evaluating quantities on cell faces. The values at cell face $j+1/2$ are shown
relative to the cell-centered values from cell $j$.

A smooth transition can of course be achieved in numerous ways, the most
popular of which is to consider a harmonic interpolation Larsen_book ;
Ramis_ComputPhysComm_2009

$\displaystyle Q_{\text{eff}}=$
$\displaystyle\,\\!\left(1/Q_{\text{Four.}}+1/Q_{\text{max}}\right)^{\\!-1}\,.$
(7)

Application of Eq. (7) reduces the effective heat flux before the transition
point (where $Q_{\text{Four.}}=Q_{\text{max}}$), which removes the blunted
shape of the thermal wave profile but also undesirably retards the propagation
of the front (Fig. 2a). Less well-known smooth schemes that do not strongly
alter the wave speed in the weakly flux-limited regime include the ‘Larson’
model

$\displaystyle Q_{\text{eff}}=$
$\displaystyle\,\\!\left(1/Q_{\text{Four.}}^{\,2}+1/Q_{\text{max}}^{\,2}\right)^{\\!-1/2}\,,$
(8)

and the ‘exponential’ model

$\displaystyle Q_{\text{eff}}=$
$\displaystyle\,Q_{\mathrm{max}}\\!\left[1-\exp\\!\left(-Q_{\mathrm{Four.}}/Q_{\mathrm{max}}\right)\right]\,.$
(9)

Both of these originate in treatments of flux-limited radiation diffusion
Mihalas_JQSRT_1982 ; Pomraning_book  and show similar behavior in how they
approach $Q_{\text{max}}$ (Fig. 2b) and accordingly lead to similar
temperature profiles for nonlinear thermal waves (Fig. 2a).

The electron and ion heat flux calculations are forced to follow Fourier’s law
(4) in Hytrac, subsequently introducing a limited form of thermal conductivity
DEIRA_code_manual ; Ramis_ComputPhysComm_2009

$\displaystyle\kappa_{\max}=$ $\displaystyle\,Q_{\max}/\\!\left|\nabla
T\right|\,.$ (10)

All the heat fluxes in the foregoing expressions (6)–(9) can therefore be
replaced with thermal conductivities, such that
$\mathbf{Q}_{\text{eff}}=-\kappa_{\text{eff}}\nabla T$ by definition and

$\displaystyle\kappa_{\text{eff}}=$
$\displaystyle\,f_{\text{intrp}}\\!\left(\kappa_{\text{tab}},\kappa_{\max}\right)\,,$
(11)

where $\kappa_{\text{tab}}$ refers to the tabulated thermal conductivity at a
given set of thermodynamic conditions and $f_{\text{interp}}$ represents one
of the interpolation functions Eqs. (6)-(9). Use of Eq. (10) crucially allows
the diffusion time step $\Delta t_{\text{dif}}=(\Delta
x)^{2}\tilde{C}_{V}/2\rho\kappa$ in the Runge-Kutta-Legendre-2 explicit super
time-stepping integration scheme Meyer_JComputPhys_2014  to be self-consistent
with the heat flow in the flux-limited regime, thereby greatly reducing run
time.

###  III.4 Evaluation of cell-face quantities

As with all codes based on the finite volume method, the heat fluxes between
EGrid cells must be calculated on the cell faces and are evaluated via
Blazek_book

$\displaystyle\left.\frac{\partial E}{\partial t}\right|_{\text{cond}}=$
$\displaystyle\,-\nabla\cdot\mathbf{Q}_{\text{eff}}\stackrel{{\scriptstyle
V\to
0}}{{=}}\frac{1}{V}\sum_{\\!\left\\{\text{faces}\right\\}}\\!\left(\mathbf{S}\cdot\mathbf{Q}_{\text{eff}}\right)_{\text{face}}\,,$
(12)

in which $V$ is the cell volume, $\mathbf{S}$ is the cell face area vector and
the summation extends over all cell faces. Thus, the effective thermal
conductivity (11) must be evaluated on the cell faces. Following Eq. (10) the
temperature gradient is approximated with a linear finite difference based on
cell-centred values, e.g. $\\!\left|\nabla
T\right|^{j+1/2}\approx\\!\left|T^{j+1}-T^{j}\right|/\delta$, where $\delta$
is the displacement between cell centres.

Assuming that the thermodynamic fields within the system are smooth and
continuous then at sufficient resolution a linear interpolation is the most
natural choice to evaluate $\kappa_{\text{tab}}$ and $Q_{\max}$ on the cell
faces. For the conduction stencil used in Hytrac this amounts to the
arithmetic mean (exemplified here using the tabulated thermal conductivity)

$\displaystyle\kappa_{\text{tab}}^{j+1/2}=$
$\displaystyle\,\frac{1}{2}\\!\left(\kappa_{\text{tab}}^{j+1}+\kappa_{\text{tab}}^{j}\right)\,,$
(13)

in which the superscripts refer to cell center (integer) and face (half-
integer) indices. Weighting the cell-centred quantities using their nonlinear
scaling on thermodynamic variables, e.g. by $T^{5/2}/Z_{i}^{*}$ for the
electron thermal conductivity in the ideal plasma regime nrl_formulary , gives
a better representation of their variation over the EGrid in comparison to Eq.
(13). For the general case one can simply self-weight the quantities,
resulting in a contraharmonic mean

$\displaystyle\kappa_{\text{tab}}^{j+1/2}=$
$\displaystyle\,\frac{(\kappa_{\text{tab}}^{j+1})^{2}+(\kappa_{\text{tab}}^{j})^{2}}{\kappa_{\text{tab}}^{j+1}+\kappa_{\text{tab}}^{j}}\,.$
(14)

For ideal states the $T^{5/2}$ scaling of the conductivity means that the heat
flux resulting from Eq. (14) is driven by the hotter of the cells in thermal
contact, i,e. the properties of the cell downstream of the temperature
gradient have no bearing on the flux it may accept.

An alternative interpolation method used in thermal transfer as a
phenomenological description of heat flow at high-thermal impedance junctions
is to take a harmonic mean

$\displaystyle\kappa_{\text{tab}}^{j+1/2}=$
$\displaystyle\,2\frac{\kappa_{\text{tab}}^{j+1}\kappa_{\text{tab}}^{0}}{\kappa_{\text{tab}}^{j+1}+\kappa_{\text{tab}}^{0}}\,.$
(15)

This model enforces continuous heat flux and temperature between cells and a
discontinuity of the temperature gradient and conductivity
Tsui_NumHeatTrans_2014  and can be thought of as complimentary to the
arithmetic mean (Fig. 2c). For the case of a solid insulator-conductor
interface Eq. (15) approximately captures the contact resistance resulting
from microscopic voids Incropera_book  and is therefore unlikely to give good
results if one of the materials is a gaseous plasma.

For the propagation of nonlinear thermal waves between regions having
conductivities differing by orders of magnitude Eq. (15) fails to recover
analytic results. For example, well-established test cases
Reinicke_PhysFluidsA_1991 ; Zeldovich_book  with zero conductivity in the
unheated part of the domain fail entirely since the resulting heat flux at the
cell boundary representing the front is exactly zero. Similar results are
shown in Fig. 2a, where the harmonic mean curve travels only $2.5\,\mathrm{\mu
m}$ from the an initial temperature discontinuity at $x_{0}=5\,\mathrm{\mu m}$
($T_{0\text{L}}=1\,\mathrm{keV}$ and $T_{0\text{R}}=293\,\mathrm{K}$) compared
to the $45-50\,\mathrm{\mu m}$ predicted using other model combinations. For
cases where a large gradient in conductivity does not exist
Shafranov_SovPhysJETP_1957 ; Kamm_LANL_2008  results using Eq. (15) perform
similarly to those using Eqs. (13) and (14).

Despite these concerns, we are motivated to include the harmonic mean option
in our study due to its inclusion in numerical solutions to the heat equation
Meyer_JComputPhys_2014  and other simulation codes used for HEDP research
vanderHolst_AstrophysJSupplSeries_2011 . Furthermore, the choice of how to
evaluate the maximum heat flux $Q_{\max}$ in Eq. (10) at cell faces is also
arbitrary and, in principle, independent to the choice applied for
$\kappa_{\text{tab}}$. In this case the same rationale behind the arithmetic
and contraharmonic mean interpolations applies, although there is no clear
analogue to the reasoning behind use of the harmonic mean approach. We
consider the latter as a means of examining flux restriction across material
interfaces under conditions where the heat flow is dominated by the free-
streaming contribution.

(a)

(b)

(c)

Figure 3:  Cell trajectory plots on the $\rho-T_{i}$ plane for the fuel
(left), coverslip (middle) and anvil (right) materials. Contours for
$\Gamma_{ii}=1$ (red), $\Gamma_{e}=1$ (blue), $D_{e}=1$ (green) and
$W=\\!\left\\{0.01,0.1,0.99\right\\}$ (orange) are shown. For the initially
solid materials (PMMA and gold) the Maxwell construction region, where the
accuracy of the FEOS model is highly uncertain, is denoted by the horizontal
isobars; There are no states in this region at any time in the simulation. The
time evolution of the simulation is denoted by the color scale (blue - early
time and dark red - late time). The initial conditions of the deuterium fuel
are shifted from true STP conditions for $P_{0}=1\,\mathrm{bar}$ due to the
lack of molecular dissociation in FEOS. Ad hoc corrections following the
linear mixing model Ross_JChemPhys_1983  have been found to make negligible
difference to the BWA profiles or final neutron yield in these simulations, as
expected for fuel initially in the gaseous phase.

###  III.5 Conduction across a material interface

With the foregoing discussion in mind one of the most important uncertainties
in our modeling is the treatment of heat flow through material interfaces;
especially for the deuterium-gold boundary, where the heated fuel is typically
several orders of magnitude more conductive and the heat flow could also be
flux-limited. If the interface is highly conductive then the internal energy
will quickly leak into the high-heat capacity anvil and quench the yield. The
same concern also holds for the later stages, when the fuel is being
compressed within the collapsing cavity. On the other hand, if the interface
is strongly insulating, which would result from use of Eq. (15), then the fuel
ions will tend to stay hot during the stagnation phase. Moreover, the
evolution of the state is coupled to electronic conduction through
equilibration, the rate of which is also influenced by the fact that the
initial nonequilibrium produced by the shock is enhanced when ion conduction
is suppressed.

Using Spitzer’s well-known expressions for the thermal conductivities we may
estimate the conditions for which the conduction across an interface between
fully-ionized deuterium (D) and some high-Z material (Z) is driven by the
conditions on one side only. Solving for the states that lead to a 10-fold
disparity, i.e. $\kappa_{e,i}^{\text{D}}\geq 10\kappa_{e,i}^{\text{Z}}$, gives

$\displaystyle T_{e\,\text{D}}\gtrsim$
$\displaystyle\,0.6772\frac{T_{e\,\text{Z}}}{Z^{*}_{\text{Z}}}\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ (\text{electrons})\,,$ (16)
$\displaystyle T_{i\,\text{D}}\gtrsim$
$\displaystyle\,2.512\,\frac{T_{i\,\text{Z}}}{\\!\left(Z_{\text{Z}}^{*}\right)^{8/5}}\,\\!\left(\frac{A_{\text{D}}}{A_{\text{Z}}}\right)^{1/5}\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ (\text{ions})\,,$ (17)

in which $A$ denotes the atomic mass and $Z_{\text{Z}}^{*}$ the effective ion
charge for transport processes in the high-Z material.

Taking $Z_{\text{Au}}^{*}\sim 6$ and $Z_{\text{PMMA}}^{*}=3$ as representative
ionization states around the time where the fusion output begins, Eqs. (16)
and (17) suggest that the electron conductivity of the fuel dominates when
$T_{e\,\text{D}}\gtrsim 0.113T_{e\,\text{Au}}$ and $T_{e\,\text{D}}\gtrsim
0.226T_{e\,\text{PMMA}}$. The equivalent criteria for the ions are
$T_{i\,\text{D}}\gtrsim 0.0818T_{i\,\text{Au}}$ and $T_{i\,\text{D}}\gtrsim
0.425T_{i\,\text{PMMA}}$. These conditions are well-fulfilled for both the
D-Au and D-PMMA interfaces in our simulations throughout the period with the
highest reactivity. A similar argument using flux-limited forms leads to
criteria which are also satisfied at all times of interest. Using Eqs.
(13)-(14) the properties of the fuel will dominate whereas using Eq. (15) the
heat flow will be dictated by the properties of the high-Z material. Thus, we
expect a strong correlation between target performance and the choice of cell
face interpolation scheme. Naturally, these conclusions hold only under the
restriction of the flux-limited model and may not properly capture the true
nature of interfacial transport.

###  III.6 Uncertainties in conduction microphysics

Aside from the uncertainties in the implementation of the conduction operator
itself one of the largest sources of uncertainty is the accuracy of the
tabulated microphysics models which underpin the energetics. In these
simulations we use the well-known Lee-More Lee_PhysFluids_1984  and Stanton-
Murillo Stanton_PhysRevE_2016  models for the conductivities of the electrons
and ions, respectively. The Lee-More model is corrected to include electron-
electron scattering following the approach of Apfelbaum
Apfelbaum_PhysPlasmas_2020 . For the electron-ion energy exchange rate we use
the $f$-sum rule approach Hazak_PhysRevE_2001 ; Gericke_JPhysConfSeries_2005 ,
which been shown to perform well compared to molecular dynamics simulations
Vorberger_HEDP_2014 . Higher-order considerations such as the coupled mode
effect Dharma-wardana_PhysRevLett_1991 ; Dharma-wardana_PhysRevE_1998 ;
Vorberger_PhysRevE_2010  can safely be ignored since the nonequilibrium states
produced in our simulations almost always have $T_{i}>T_{e}$. All the
microphysics models are driven by the ionization predicted by FEOS for
consistency with the EoS Faik_ComputPhysCommun_2018 .

Whilst each of these models include corrections for electron degeneracy and
strong ion coupling they are still not expected to be very accurate under
conditions in warm dense matter (WDM) regime. Of particular concern is the
electron thermal conductivity as this is strongly influenced by many-body
effects such as screening, structural effects and the ionization equilibrium
of high-Z systems Ichimaru_PhysRevA_1985b ; Kitamura_PhysRevE_1995 ;
Hu_PhysRevE_2014b ; Hu_PhysPlasmas_2016 ; Starrett_HEDP_2020 ;
Shaffer_PhysRevE_2020 .

A simple but informative a priori indicator for model accuracy is encapsulated
by the WDM parameter Murillo_PhysRevE_2010

$\displaystyle W(\rho,T_{e})=$
$\displaystyle\,S(\Theta)S(\Gamma_{ee}),\leavevmode\nobreak\
\leavevmode\nobreak\ S(x)=2/(x+1/x)\,,$ (18)

which provides a simple measure of the nonideality of thermodynamic states
from the perspective of theoretical modeling. In Eq. (18), the usual
definitions of the degeneracy parameter, $\Theta$, and electron coupling
parameter, $\Gamma_{ee}$, are used Murillo_PhysRevE_2010 . Cell trajectory
plots extracted from the reference simulation (Fig. 3) show that the deuterium
fuel is mostly ensconced within the ideal plasma regime, where $D_{e}\ll 1$
and $\Gamma_{ee,ii}\ll 1$, such that $W\ll 1$. Although at late time, where
the neutron emission rate is highest, the degree of nonideality in the fuel is
seen to steadily increase. Significantly larger uncertainties can be expected
for the properties of the coverslip and anvil materials as they produce states
with nonideal electrons, $W\sim 1$, and strongly coupled ions
$\Gamma_{ii}\gtrsim 1$.

We account for the potential uncertainties in the microphysics by assessing
the impact of a range of scaling factors, $s$, covering two orders of
magnitude; $s=\\!\left\\{0.1,0.5,2,10\right\\}$. These may be applied
uniformly across phase space, referred to as a ‘blanket approach’, or
conversely by a ‘targeted approach’ which utilizes the WDM parameter Eq. (18)
to focus their application toward regions with high uncertainty. The effective
scaling parameter for a particular point in $\rho-T_{e}$ space is then

$\displaystyle s^{\prime}(\rho,T_{e})=$
$\displaystyle\,1+\\!\left(s-1\right)W\\!\left(\rho,T_{e}\right)\,.$ (19)

The targeted form of the scaling factor crucially ensures that the tabulated
microphysics are not unreasonably distorted under conditions where uncertainty
in their forms is believed to be negligible, e.g. in the high-temperature,
low-density limit where the well-established Spitzer-type models apply.

###  III.7 Kinetic (nonlocal) effects

A simple estimate of the degree of locality with regard to collision-driven
transport processes is the Knudsen number

$\displaystyle N_{\text{K}j}=$ $\displaystyle\,\frac{\lambda_{j}}{L}\,.$ (20)

Here the collisional mean free path $\lambda_{j}$ of species $j=e,i$ can be
estimated from simple expressions spitzer_book  throughout the simulation
domain as a function of local thermodynamic variables. For the scale length,
$L$, at a particular location in the system, we take the harmonic mean of
inverse logarithmic derivatives

$\displaystyle L\approx$
$\displaystyle\,N_{X}\\!\left(\sum_{X}\frac{\\!\left|\nabla
X\right|}{X}\right)^{-1}\,,$ (21)

where $X=\\!\left\\{T_{i},T_{e},\rho,Z_{i}^{*}\right\\}$ represents the set of
variables which most strongly influence the collisional mean free path and
$N_{X}=4$ in this case. This definition picks up both the change in density at
a material interface and the strong gradient in heating at shock fronts and is
similar to the scheme discussed by Taitano et al. Taitano_PhysPlasmas_2018 .

(a)

(b)

(c)

Figure 4:  Heat maps of ion temperature (left-hand side) and ion Knudsen
number (right-hand side) computed using Eqs. (20) and (21) for the reference
simulation discussed in Section II.2 for different times throughout the cavity
implosion: Just prior to stagnation of the incident shock at the tip of the
cavity (left); After stagnation, where the reflection from the tip drives up
into the collapsing cavity (middle); During the quasi-isothermal compression
phase (right). The Knudsen number is initially of order unity close to the
hottest region but quickly falls to a level at which local ion conduction is
appropriate. In these plots the bright green contour denotes the Lagrangian
grid representing the material interface.

For the times relevant to the fusion output from our simulations
$N_{\text{K}e}\ll 1$, such that electronic conduction is likely to be well-
described by purely local diffusion (4). On the other hand, we find
$N_{\text{K}i}\sim 0.1-1$ in a well-resolved layer within the fuel in
proximity to the D-Au interface, suggestive of strongly nonlocal ionic
transport into the anvil material at the point of shock stagnation where the
highest ion temperatures occur (Fig. 4). This nonlocal layer dissipates
quickly after shock reflects from the cavity tip.

Since the conditions in the fuel throughout the periods where nonlocal effects
may be significant are conducive to a large contribution from the ions we
expect this to have a strong influence on the fuel energetics and, thus, the
fusion performance. Moreover, the reduction of the thermal reactivity
Molvig_PhysRevLett_2012 ; Albright_PhysPlasmas_2013 ; Kagan_PhysRevLett_2015
is likely to be substantial for the initial shock stagnation event. This is
presently being investigated and will reported in a future publication.

##  IV Results and discussion

To assess the impact of the configuration changes relative to the reference
simulation (Tab. 1) we track three simple metrics which are sensitive to both
the dynamics of the incident shock and the fuel energetics throughout the
period of neutron emission:

  1. 1.

Time of first reflection of the incident shock from the axis of the
simulation, $t_{\text{axis}}$

  2. 2.

Maximum BWA ion temperature during the stagnation phase, $T_{\max}$.

  3. 3.

Mean value of BWA ion temperature during the compression phase,
$T_{\text{av}}$.

These are presented as as percentage changes relative to the corresponding
values obtained from the reference simulation, i.e.
$\delta_{i}=100\times\\!\left(M_{i}/M_{i}^{\text{ref}}-1\right)$, where
$M_{i}$ stands for the relevant metric ($i=\\{1,2,3\\}$). We originally
considered several other metrics such as the full-width-half-maximum duration
of the stagnation phase and the standard deviation of the BWA ion temperature
in the compression phase; both of these proved unreliable in producing stable
trends over the wide range of perturbations to the reference configuration.
Most importantly, we have chosen not to track the total neutron yield
resulting from each configuration change. This is principally because it
scales so nonlinearly with ion temperature that any transient artefacts
resulting from small pockets of fuel caught in the complex structures formed
lead to spurious results that are both essentially impossible to correct for
and distort otherwise meaningful and insightful trends. This is an unfortunate
but inescapable consequence of the very low yields ($Y_{\text{n}}\gtrsim 100$)
predicted for this target.

Variable | Value
---|---
Electron flux limiter interp. | Exponential
Ion flux limiter interp. | Exponential
Electron flux limiter coef. | 0.05
Ion flux limiter coef. | 0.5
Conductivity cell face interp. | Arithmetic mean
Max heat flux cell face interp. | Arithmetic mean
Table 1:  Configuration options of Hytrac’s flux-limited conduction operator
used in the reference simulation.  Figure 5:  Effect of scaling factors
applied on the conduction microphysics models featured the reference
simulation; electron thermal conductivity (left-hand column), ion thermal
conductivity (middle column) and electron-ion energy exchange rate (right-hand
column). These are varied in the deuterium fuel (top row), gold anvil (middle
row) and PMMA coverslip (bottom row). The three metrics described in the text
are plotted in each panel as a percentage change relative to the reference
case: Metric 1 (blue markers/line); Metric 2 (orange markers/line); Metric 3
(green marker/line). The horizontal shaded band in each panel between
$\delta=\pm 10\%$ represents the maximum impact relative to the reference
simulation that we (arbitrarily) deem significant. In each panel the blanket
(filled marker connected by solid lines) and targeted (open markers connected
by dashed lines) approaches to applying the scaling factors are used.

We first examine the impact of the scaling factors applied to the components
of the conduction microphysics; the electron and ion thermal conductivities
and the electron-ion energy exchange rate. All three materials which directly
contribute to the energetics of the collapsing cavity (the deuterium fuel, the
PMMA coverslip and the gold anvil) are considered. For each material, scaling
factors in the range $s=\\!\left\\{0.1,0.5,2,10\right\\}$ have been applied
(individually) using the ‘blanket’ and ‘targeted’ manner as described in
Section III.6. To simplify the analysis we consider only order-one
perturbations, i.e. only a single aspect of the system is perturbed at a time.
This has the disadvantage that interactions between uncertainties are ignored,
which is surely important in a two-temperature, conduction-dominated plasma.
On the other hand, establishing the basic trends is of prime importance and
lays the foundations for more complex, coupled sensitivity studies in future.

Strong trends are found with respect to all the fuel microphysics using the
blanket approach to applying the scaling factors (filled markers connected by
solid lines in Fig. 5), particularly in the maximum BWA ion temperature during
stagnation. This verifies the assertion that transport and equilibration play
a crucial role in the energetics of the fuel at stagnation. However, the
strongest correlation is between the electronic conductivity of the fuel and
the average temperature in the compression phase, with substantially less
sensitivity to the ionic conductivity and equilibration rate. Electron
conduction is therefore the most important loss mechanism during the period
over which the majority of the neutron output is expected. This can be readily
understood since the plasma is essentially in thermal equilibrium during this
phase, such that the relative importance of conduction is dictated by the
ratio $\kappa_{e}/\kappa_{i}\sim(m_{i}/m_{e})^{1/2}\gg 1$. Weaker correlations
are evident for the properties of the anvil and coverslip. The lack of any
clear influence of any of the scaling factors on the arrival time of the first
shock on the simulation axis (Metric 1) can be taken to mean that the early
time dynamics of this target are mostly hydrodynamic in nature.

In all cases the sensitivity of the metrics to any of the scaling factors
drops when the targeted approach (19) is used (open markers connected by
dashed lines in Fig. 5). In this more realistic case none of the perturbations
lead to a change beyond $\pm 10\%$, which we arbitrarily define to constitute
a ‘significant’ deviation from the reference case. For the fuel this is
unsurprising given the largely ideal states produced (Fig. 3a). Although
larger values of the scaling factor can be substantiated for the initially
solid materials the results show that their thermophysical properties do not
greatly influence target performance. The only exception is the possibility
that the electron thermal conductivity of gold is underestimated by the LM
model by a factor of 10 in the WDM regime. Although we have not been able to
find any data in the literature to substantiate such a large discrepancy it
could easily be accommodated by the general paucity of the model under WDM
conditions. Firmer conclusions will only be able to be reached when more
accurate thermal conductivity tables are available.

Figure 6:  Same as in Fig. 5, but considering the influence of the choice of
interpolation scheme for evaluating the tabulated conductivities (left) and
max heat fluxes (right) on cell faces using cell-centered quantities (top
row), the electron (left) and ion (right) flux limiter coefficients (middle
row) and the interpolation method which applies the limited form of the
electron (left) and ion (right) heat flux (bottom row).  Figure 7:  Comparison
of the BWA ion temperature (top panel) and cumulative yield (bottom panel)
from the reference case to those from simulations with the three most
significant changes to the stagnation temperature. Note the change in the
temporal scale on the $x$-axis between the top and bottom panels.

Another factor which explains the lack of observed sensitivity to the
properties of the anvil and coverslip is the dominant role of the fuel in
determining the heat loss through the material interfaces, which can be
understood since the conditions encapsulated by Eqs. (16) and (17) are well-
fulfilled and the interfacial flux in these simulations is constructed with
the arithmetic mean approach (13). Larger uncertainties consistent with the
nonideal conditions in the anvil might therefore be expected if the interface
is better described by a harmonic mean approximation. This has not been
repeated here since the indications of the test cases discussed in Sec. III.4
are that using Eq. (15) will not properly capture nonlinear thermal wave
propagation through large conductivity gradients. Nevertheless, it is still
interesting to note its wider impacts on the target performance and in the
absence of direct verification or validation evidence it cannot presently be
completely ruled as a modeling option.

As shown in the top row of Fig. 6, where the interpolation of cell-centered
quantities to the cell faces is performed using contraharmonic and harmonic
averages, Eqs. (14) and (15), the impact of using the harmonic mean in
particular is indeed significant. For the harmonic mean case the insulation of
the fuel resulting from the dominance of the non-conductive anvil in setting
the conduction flux increases the peak stagnation temperature by a factor of
two to three ($T_{\max}\gtrsim 1-1.3\,\mathrm{keV}$) compared to the reference
case. The result is a single ‘flash’ of neutrons, which constitutes almost all
of the yield from the target as the contribution from the state of the plasma
created in the compression phase is essentially unaffected and does not reach
a competitive fusion rate (Fig. 7).

This fundamental change in the character of the neutron production is observed
if either the unlimited (i.e. the tabulated data) or limited (10) components
of the effective thermal conductivities are changed independently, with
roughly $\sim 30\%$ higher ion temperatures achieved when only the maximum
heat flux is constructed on the cell face using Eq. (15). The corresponding
amplification in the yield is much larger, being roughy $\sim 6$ times larger
than for the case where only the tabulated conductivity is harmonically
averaged. Taken together, these results are indicative of the ionic heat flow
being moderately flux-limited, such that both contributions make substantial
contributions to the total effective flux. Usage of the harmonic mean also
leads to a faster initial shock, as shown by the earlier arrival time at the
axis, and increases the temperature of the state formed by reflection from the
cavity wall just prior to stagnation. The average temperature during the
compression phase is, however, very weakly modified since at this stage most
of the conduction occurs within the fuel where temperature (conductivity)
gradients are shallow. In contrast, using the contraharmonic mean for either
component of the effective conductivity has a negligible effect on any of the
metrics.

The target operation is found to be similarly affected without using the
harmonic mean by suppressing ionic conduction through the flux limiter
coefficient and to a far lesser degree by the choice of interpolation function
used to apply the limited form (see right-hand panels of the bottom and middle
row of Fig. 6, respectively). Setting $\alpha_{i}=\alpha_{e}=0.05$, which
amounts to a ten-fold reduction in the maximum ion heat flux from the
reference simulations, again leads to a peak BWA ion temperature nearly twice
as hot as the reference case during stagnation and similarly increased neutron
total yield from stagnation, whilst the other metrics are essentially
unaffected, but is substantially reduced relative to the purely diffusive
expression (4). This is consistent with the findings from the scaling factor
study where the role of ion conduction during the compression phase was shown
to be negligible. The fact that there is very little change to any of the
metrics if the fraction of the maximum ion heat flux is allowed to increase
arbitrarily, i.e. for $\alpha_{i}\geq 1$, again suggests that the ionic heat
flow is not strongly flux-limited. This is further backed up by the
observation that switching to the softer harmonic interpolation method (7)
acts to slightly increase the stagnation temperature by cutting the effective
heat flux, whereas the min-capped approach (6) slightly enhances it. The ratio
$Q_{\text{eff}}/Q_{\max}$ is then close to unity (see Fig. 2b) since in the
strongly flux-limited regime one would expect a much larger effect, especially
from using harmonic interpolation. In contrast, the evidence does not suggest
that electron heat flow is flux-limited as no significant sensitivities are
observed for any of the configuration changes considered, in agreement with
the expectations noted in Sec. III.7.

Despite the indication being that the ionic heat flux is quasi-local it is
important to note that only separate perturbations to the operator
configuration have been examined, meaning that the potential for large changes
to the performance caused by coupling of multiple configuration changes must
be assessed before stronger conclusions can be drawn. Taking into account the
observed strong sensitivities to the treatment of conduction at the interface
and the lack of any compelling evidence for the reasonableness of a single,
universal value of $\alpha_{i}\gg 0.05$, as well as the fact that high ion
Knudsen numbers of order $0.1-1$ are predicted around stagnation (Fig. 4), it
is difficult to conclude that the reference case is suitably predictive,
especially from the perspective of being able to quantitatively compare
results for the total yield with the reported experimental values from Ref.
Derentowicz_JourTechPhys_1977 . Nevertheless, we note that none of the changes
considered here result in a total yield close to the lowest reported yields of
Derentowicz et al., with the highest total yield being roughly fives smaller.
The work required to bring together all these threads in a properly coupled
manner is presently ongoing, as are efforts to develop our understanding of
the possible role and impact of kinetic ion transport on these targets and the
influence of the other radiation-hydrodynamic phenomena left out of this work.

##  V Conclusions

We have undertaken an initial examination of the influence of flux-limited
thermal conduction modeling on the robustness and performance of the
uniaxially-driven conical fusion targets used in the experiments of
Derentowicz et al. Derentowicz_JourTechPhys_1977 . This kind of target is of
interest to FLF as it is one of the few examples of fusion targets in the
literature with observable neutron yield, where the initial collapse is driven
by a single planar shock produced by conditions within reach of our current
experimental facilities. The attention given to thermal conduction in this
paper is justified since other phenomena such as radiation transport and
viscous effects are not expected to play as important a role in determining
the energetics of the fusion fuel.

The uncertainties which principally affect conduction loss are associated with
the related plasma microphysics and the degree to which the transport is
nonlocal. To assess the impact of uncertainties in the microphysics we have
undertaken scaling factor studies in which the tabulated models are scaled
both uniformly, i.e. a blanket approach, and with a more realistic method
using the WDM parameter Murillo_PhysRevE_2010  to target regions of density-
temperature space with larger theoretical uncertainty. The reference case was
produced using an idealized model that was rigorously converged on the basis
of root-mean-square differences between time-dependent profiles of burn-
weighted average (BWA) plasma parameters from successively higher resolution
simulations. Simple metrics which capture the dynamics of the initial stages
of the cavity collapse and the energetics of the fusion fuel over the time of
neutron emission were defined, from which direct comparison can be made
between different simulations and the reference case. Uniform application of
the scaling factors showed that electron conduction in the fuel during the
compression phase most strongly affects the target performance. With the more
realistic targeted approach we found very little influence on any of the
metrics. This is due to the fact that it is the properties of the fuel that
determine its energetics in the compression phase, which remains strongly
ideal, and therefore minimally uncertain, throughout the simulations. We found
that the role of the nonideal states expected to be produced in the plastic
coverslip and gold anvil is minimised due to the numerical treatment of the
heat flux at material interfaces.

Considering alternative interpolation schemes for evaluating cell-centered
quantities at cell faces leads to significant changes to the predicted
conditions and also in the character of the neutron emission. In particular,
we found that the largest influence on target performance comes from changing
the cell-face averaging from the arithmetic mean used in the reference case to
a harmonic mean. Although the latter is not thought to be an especially
physical choice as it was shown to be incapable of propagating nonlinear
thermal waves between regions of substantially different conductivities, it
cannot be explicitly ruled out without direct validation. With this change the
suppression of conduction losses substantially impacts all of the metrics,
crucially changing the fusion output from one dominated by heating due to
quasi-isentropic compression during the cavity collapse to one in which a
single neutron flash is produced by the reflection of the incident shock from
the cavity tip during stagnation. A similar change was found to result from
restricting the ionic heat flux in the reference case from the stagnated fuel
to the anvil through a reduced flux limiter (set equal to the widely used
canonical value used in flux-limited electron conduction Atzeni_book ;
Bell_PhysRevLett_1981 ). Based on the results we conclude that the ionic heat
flow is likely to be moderately flux-limited, i.e. nonlocal kinetic effects
may be important, whereas the electronic heat flow is evidently more
diffusive.

If the ionic conduction is indeed minimally restricted, either because
Hytrac’s treatment of the heat flux across material interfaces is reasonable
or because a small value of the ionic flux limiter coefficient is
unreasonable, then our reference case simulations should be relatively robust.
Consequently, the question of why the predicted total yield of these targets
is several orders of magnitude less than those reported in Ref.
Derentowicz_JourTechPhys_1977  remains an open one. Moreover, we noted that
none of the singular changes to the reference case considered resulted in
neutron yield amplified by more than an order of magnitude. More extensive
sensitivity studies incorporating multiple, simultaneous, interacting
configuration changes as well as more elements of Hytrac’s multi-physics
model, e.g. radiation transport and viscous effects, are needed to develop a
clearer picture.

Finally, we note that the findings related to the potentially important role
of ionic conduction are interesting as they suggest that targets similar to
those studied here may provide an experimental platform for accessing states
of matter in the quasi-local regime, in which energy transport is close to the
transition between local diffusion and the kinetic transport regime. With
suitable diagnostics to provide time-resolved histories of neutron production
and electron and ion temperatures from the stagnated plasma, it may be
possible to experimentally distinguish between predictions based on simple ion
transport models (such as flux-limited conduction) and state-of-the-art
kinetic models, whilst minimizing uncertainty associated with electron
transport and equilibration. Such data may prove to be invaluable for
benchmarking integrated simulation codes for ICF.

## Data availability

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

## Acknowledgements

The authors gratefully acknowledge R. King and J. Gardiner for providing HPC
support and A. Venskus and J. Herring for invaluable contributions to the
robustness and stability of Hytrac. We also recognize insightful conversations
with A. Crilly and G. Kagan (Imperial College London) regarding equation of
state and ionic transport modeling.

## References

  * [1] H. Derentowicz, S. Kaliski, and Z. Ziółkowski.  Generation of fusion neutrons in a deuterium filled cone by means of explosive implosion of polythene shell. Part 1. Theoretical estimations.  Journal of Technical Physics, 18:465–471, 1977.
  * [2] J. A. Gaffney, D. Clark, V. Sonnad, and S. B. Libby.  Development of a Bayesian method for the analysis of inertial confinement fusion experiments on the NIF.  Nuclear Fusion, 53:073032, 2013.
  * [3] J. A. Gaffney, D. Clark, V. Sonnad, and S. B. Libby.  Bayesian inference of inaccurancies in radiation transport physics from inertial confinement fusion experiments.  High Energy Density Physics, 9:457, 2013.
  * [4] J. Melvin, H. Lim, V. Rana, B. Cheng, J. Glimm, D. H. Sharp, and D. C. Wilson.  Sensitivity of inertial confinement fusion hot spot properties to the deuterium-tritium fuel adiabat.  Physics of Plasmas, 22:022708, 2015.
  * [5] T. M. Guymer, A. S. Moore, J. Morton, J. L. Kline, S. Allan, N. Bazin, J. Benstead, C. Bentley, A. J. Comley, J. Cowan, K. Flippo, W. Garbett, C. Hamilton, N. E. Lanier, K. Mussack, K. Obrey, L. Reed, D. W. Schmidt, R. M. Stevenson, J. M. Taccetti, and J. Workman.  Quantifying equation-of-state and opacity errors using integrated supersonic diffusive radiation flow experiments on the National Ignition Facility.  Physics of Plasmas, 22:043303, 2015.
  * [6] I. K. Krasyuk, A. Yu. Semenov, and A. A. Charakhch’yan.  Use of conic targets in inertial confinement fusion.  Quantum Electronics, 35(9):769, 2005.
  * [7] A. Shyam and M. A. H. A. D. E. V. A. Srinivasan.  Investigation of the feasibility of conical target compression experiments using exploding foil driven hyper-velocity liners.  Atomkernenergie Kerntechnik, 44:196, 1984.
  * [8] V. I. Vovchenko, A. S. Goncharov, Iu. S. Kasianov, O. V. Kozlov, I. K. Krasiuk, A. A. Maliutin, M. G. Pastukhov, P. P. Pashinin, and A. M. Prokhorov.  Generation of thermonuclear neutrons by laser action on a conical target.  Journal of Experimental and Theoretical Physics Letters, 26:628, 1977.
  * [9] R. J. Mason, R. J. Fries, and E. H. Farnum.  Conical targets for implosion studies with a $\mathrm{CO}_{2}$ laser.  Applied Physics letters, 34:14, 1979.
  * [10] S. L. Bogolyubskii, B. P. Gerasimov, V. I. Liksonov, A. P. Mikhailov, Y. P.Popov, L. I. Rudakov, A. A. Samarskii, and V. P. Smirnov.  Thermonuclear neutron yield from a plasma compressed by a shell.  Journal of Experimental and Theoretical Physics Letters, 24:182, 1976.
  * [11] Jon Larsen.  Foundations of High-Energy-Density Physics.  Cambridge University Press, 2017.
  * [12] F. Graziani, M. P. Dejarlais, R. Redmer, and S. B. Trickey.  Frontiers and Challenges in Warm Dense Matter: Lecture Notes in Computational Science and Engineering.  Springer International Publishing Switzerland, 2014.
  * [13] O. Larroche.  Kinetic simulations of fuel ion transport in ICF target implosions.  European Physics Journal D, 27:131, 2003.
  * [14] O. Larroche.  Ion Fokker-Planck simulation of $\mathrm{D}-{}^{3}\mathrm{He}$ gas target implosions.  Physics of Plasmas, 19:122706, 2012.
  * [15] W. T. Taitano, L. Chacón, and A. N. Simakov.  An adaptive, implicit, conservative, $1D-2V$ multi-species Vlasov-Fokker-Planck multi-scale solver in planar geometry.  Journal of Computational Physics, 365:173, 2018.
  * [16] H. G. Rinderknecht, P. A. Amendt, S. C. Wilks, and G. Collins.  Kinetic physics in ICF: presentunderstanding and future directions.  Plasma Physics and Controlled Fusion, 60:064001, 2018.
  * [17] G. P. Schurtz, Ph. D. Nicolaï, and M. Busquet.  A nonlocal electron conduction model for multidimensional radiation hydrodynamics codes.  Physics of Plasmas, 7:4238, 2000.
  * [18] D. Cao, G. Moses, and J. Delettrez.  Improved non-local electron thermal transport model for two-dimensional radiation hydrodynamics.  Physics of Plasmas, 22:082308, 2015.
  * [19] N. M. Hoffman, G. B. Zimmerman, K. Molvig, H. G. Rinderknecht, M. J. Rosenberg, B. J. Albright, A. N. Simakov, H. Sio, A. B. Zylstra, M. Gatu Johnson, F. H. Séguin, J. A. Frenje, C. K. Li. R. D. Petrasso, D. M. Higdon, G. Srinivasan, V. Yu. Glebov, C. Stoeckl, W. Seka, and T. C. Sangster.  Approximate model for the ion-knetic regime in inertial-confinement-fusion capsule implosions.  Phyics of Plasmas, 22:052707, 2015.
  * [20] J. P. Brodrick, R. J. Kingham, M. M. Marinak, M. V Patel, A. V. Chankin, J. T. Omotani M. V. Umansky, D. Del Sorbo, B. Dudson, J. T. Parker, G. D. Kerbel, M. Sherlock, and C. P. Ridgers.  Testing nonlocal models of electron thermal conduction for magnetic and inertial confinement fusion applications.  Physics of Plasmas, 24:092309, 2017.
  * [21] M. Sherlock, J. P. Brodrick, and C. P. Ridgers.  A comparison of non-local electron transport model for laser-plasmas relevant to inertial confinement fusion.  Physics of Plasmas, 24:082706, 2017.
  * [22] R. J. Henchen, M. Sherlock, W. Rozmus, J. Katz, D. Cao, J. P. Palastro, and D. H. Froula.  Observation of Nonlocal Heat Flux Using Thomson Scattering.  Physical Review Letters, 121:125001, 2018.
  * [23] R. J. Henchen, M. Sherlock, W. Rozmus, J. Katz, P. E. Masson-Laborde, D. Cao, J. P. Palastro, and D. H. Froula.  Measuring heat flux from collective Thomson scattering with non-Maxwellian distribution functions.  Physics of Plasmas, 26:032104, 2019.
  * [24] A. R. Bell, R. G. Evans, and D. J. Nicholas.  Electron energy transport in steep temperature gradients in laser-produced plasmas.  Physical Review Letters, 46:243, 1981.
  * [25] N. B. Meezan, D. T. Woods, N. Izumi, H. Chen, H. A. Scott, M. B. Schneider, D. A. Liedahl, O. S. Jones, G. B. Zimmerman, J. D. Moody, O. L. Landen, and W. W. Hsing.  Evidence of restricted heat transport in National Ignition Facility Hohlraums.  Physics of Plasmas, 27:102704, 2020.
  * [26] D. Jian, B. Fix, J. Glimm, J. Xicheng, L. Xiolin, L. Yuanhua, and W. Lingling.  A simple package for front tracking.  Journal of Computational Physics, 213(2):613, 2006.
  * [27] M. Gittings, R. Weaver, M. Clover, T. Betlach, N. Byrne, R. Coker, E. Dendy, R. Hueckstaedt, K. New, W. R. Oakes, D. Ranta, and R. Stefan.  The RAGE radiation-hydrodynamic code.  Computational Science and Discovery, 1:015005, 2008.
  * [28] E. L. Vold, R. M. Rauenzahn, C. H. Aldrich, K. Molvig, A. N. Simikov, and B. M. Haines.  Plasma transport in an Eulerian AMR code.  Physics of Plasmas, 24:042702, 2017.
  * [29] D. Mihalas and R. Weaver.  Time-dependent radiative transfer with automatic flux limiting.  Journal of Quantitative Spectroscopy and Radiative Transfer, 28(3):213, 1982.
  * [30] Kristopher McGlinchey.  Radiation-Hydrodynamics Simulations of the Impact of Instabilities and Asymmetries on Inertial Confinement Fusion.  PhD thesis, Imperial College London, 2017.
  * [31] S. Faik, A. Tauschwitz, and I. Iosilevskiy.  The equation of state package FEOS for high energy density matter.  Computer Physics Communications, 227:117, 2018.
  * [32] H.-S. Bosch and G. M. Hale.  Improved fomrulas for fusion cross-sections and thermal reactivities.  Nuclear Fusion, 32(4):611, 1992.
  * [33] S. Atzeni and J. Meyer-ter-Vehn.  The Physics of Inertial Fusion.  Oxford University Press, 2004.
  * [34] C. R. D. Brown, D. J. Hoarty, S. F. James, D. Swatton, S. J. Hughes, J. W. Morton, T. M. Guymer, M. P. Hill, D. A. Chapman, J. E. Andrew, A. J. Comley, R. Shepherd, J. Dunn, H. Chen, M. Schneider, G. Brown, P. Beiersdorfer, and J. Emig.  Measurements of electron transport in foils irradiated with a picosecond time scale laser pulse.  Physical Review Letters, 106:185003, 2011.
  * [35] D. J. Hoarty, P. Allan, S. F. James, C. R. D. Brown, L. M. R. Hobbs, M. P. Hill, J. W. O. Harris J. Morton, M. G. Brookes, R. Shepherd, J. Dunn, H. Chen, E. Von Marley, P. Beiersdorfer, H. K. Chung, R. W. Lee, G. Brown, and J. Emig.  Observations of the effect of ionization-potential depression in hot dense plasma.  Physical Review Letters, 110:265003, 2013.
  * [36] C.-K. Huang, K. Molvig, B. J. Albright, E. S. Dodd, E. L. Vold, G. Kagan, and N. M. Hoffman.  Study of the ion kinetic effects in ICF runaway burn using a quasi-1D hybrid model.  Physics of Plasmas, 24:022704, 2017.
  * [37] M. Rosen, H. Scott, D. Hinkel, E.Williams, D. Callahan, R. Town, L. Divol, P. Michel, W. Kruer, L. Suter, R. London, J. Harte, and G. Zimmerman.  The role of a detailed configuration accounting (DCA) atomic physics package in explaining the energy balance in ignition-scale hohlraums.  High Energy Density Physics, 7:180, 2011.
  * [38] J. T. Larsen and S. M. Lane.  HYADES - A plasma hydrodynamics code for dense plasmas.  Journal of Quantitative Spectroscopy and Radiative Transfer, 51:179, 1994.
  * [39] M. Basko.  DEIRA: A 1-D 3-T hydrodynamic code for simulating ICF targets driven by fast ion beams.  Technical report, Institute for Theoretical and Experimental Physics, Moscow.
  * [40] J. D. Huba.  NRL Plasma Formulary.  The Office of Naval Research, 2016.
  * [41] R. Ramis, J. Meyer-ter-Vehn, and J. Ramírez.  MULTI2D – a computer code for two-dimensional radiation hydrodynamics.  Computer Physics Communications, 180:977, 2009.
  * [42] G. C. Pomraning.  The Equations of Radiation Hydrodynamics.  Courier Corporation, 2005.
  * [43] C. D. Meyer, D. S. Balsara, and T. Aslam.  A stabilized Runge-Kutta-Legendre method for explicit super-time-stepping of parabolic and mixed equations.  Journal or Computational Physics, 257:594, 2014.
  * [44] J. Blazek.  Computational Fluid Dynamics: Principles and Applications.  Butterworth-Heinemann, 2005.
  * [45] Tsui Y.-Y., Lin S.-W., and Ding K.-J.  Modeling of heat transfer across the interface in two-fluid flows.  Numerical Heat Transfer, Part B: Fundamentals, 66:162, 2014.
  * [46] F. P. Incropera and D. P. De Witt.  Fundamentals of Heat and Mass Transfer.  Wiley, 2002.
  * [47] P. Reinicke and J. Meyer-ter-Vehn.  The point explosion with heat conduction.  Physics of Fluids A: Fluid Dynamics, 3:1807, 1991.
  * [48] Y. B. Zeldovich and Y. P. Raizer.  Physics of Shock Waves and High Temperature Hydrodynamic Phenomena.  Dover, Mineola, NY, 2002.
  * [49] V. D. Shafranov.  The structure of shock waves in a plasma.  Soviet Physics Journal of Experimental and Theoretical Physics, 5:1183, 1957.
  * [50] J. R. Kamm, J. S. Brock, S. T. Brandon, D. L. Cotrell, B. Johnson, P. Knupp, W. J. Rider, T. G. Truncano, and V. G. Weirs.  Enhanced verifiation testsuite for physics simulation code.  Technical Report LA-14379, Los Alamos National Laboratory, 2008.
  * [51] B. van der Holst, G. Tóth, I. V. Sokolov, K. G. Powell, J. P. Holloway, E. S. Myra, Q. Stout, M. L. Adams, J. E. Morel, S. Karni, B. Fryxell, and R. P. Drake.  CRASH: A block-adaptive-mesh code for radiative shock hyrdodynamics \- Implementation and verification.  The Astrophysical Journal Supplement Series, 23:194, 2011.
  * [52] M. Ross, F. H. Ree, and D. A. Young.  The equation of state of molecular hydrogen at very high density.  The Journal of Chemical Physics, 79:1487, 1983.
  * [53] Y. T. Lee and R. M. More.  An electron conductivity model for dense plasmas.  Physics of Fluids, 27:1273, 1984.
  * [54] L. G. Stanton and M. S. Murillo.  Ionic transport in high-energy-density matter.  Physical Review E, 93:043203, 2016.
  * [55] E. M. Apfelbaum.  The calculations of thermophysical properties of low-temperature gallium plasma.  Physics of Plasmas, 22:092703, 2020.
  * [56] G. Hazak, Z. Zinamon, Y. Rosenfeld, and M. W. C. Dharma-wardana.  Temperature relaxation in two-temperature states of dense electron-ion systems.  Physical Review E, 64:066411, 2001.
  * [57] D. O. Gericke.  Kinetic approach to temperature relaxation in dense plasmas.  Journal of Physics: Conference Series, 11:111, 2005.
  * [58] J. Vorberger and D. O. Gericke.  Comparison of electron-ion energy transfer in dense plasmas obtained from numerical simulations and quantum kinetic theory.  High Energy Density Physics, 10:1, 2014.
  * [59] M. W. C. Dharma-wardana.  Nature of coupled-mode contributions to hot electron relaxation in semiconductors.  Physical Review Letters, 66:197, 1991.
  * [60] M. W. C. Dharma-wardana and F. Perrot.  Energy relaxation and the quasiequation of state of a dense two-temperature nonequilibrium plasma.  Physical Review E, 3705:58, 1998.
  * [61] J. Vorberger, D. O. Gericke, Th. Bornath, and M. Schlanges.  Energy relaxation in dense, strongly coupled two-temperature plasma.  Physical Review E, 81:046404, 2010.
  * [62] S. Ichimaru and S. Tanaka.  Theory of interparticle correlations in dense, high-temperature plasmas. V. Electric and thermal conductivities.  Physical Review A, 32:1790, 1985.
  * [63] H. Kitamura and S. Ichimaru.  Electric and thermal resistivities in high-Z plasmas.  Physical Review E, 51:6004, 1995.
  * [64] S. X. Hu, L. A. Collins, V. N. Goncharov, T. R. Boehly, R. Epstein, R. L. McCrory, and S. Skupsky.  First-principles opacity table of warm dense deuterium for inertial-confinement-fusion applications.  Physical Review E, 90:033111, 2014.
  * [65] S. X. Hu, L. A. Collins, V. N. Gocharov, J. D. Kress, R. L. McCrory, and S. Skupsky.  First-principles investigations on ionization and thermal conductivity of polystyrene for inertial confinement fusion applications.  Physics of Plasmas, 23:042704, 2016.
  * [66] C. E. Starrett, N. R. Shaffer, D. Saumon, R. Perriot, T. Nelson, L. A. Collins, and C. Ticknor.  Model for the electrical conductivity in dense plasma mixtures.  High Energy Density Physics, 36:100752, 2020.
  * [67] N. R. Shaffer and C. A. Starrett.  Model of electron transport in dense plasmas spanning temperature regimes.  Physical Review E, 101:053204, 2020.
  * [68] M. S. Murillo.  X-ray Thomson scattering in warm dense matter at low frequencies.  Physical Review E, 81:036403, 2010.
  * [69] L. Spitzer.  Physics of Fully Ionized Gases.  Interscience, New York, NY, 1962.
  * [70] W. T. Taitano, A. N. Simakov, L. Chacón, and B. Keenan.  Yield degradation in inertial-confinement-fusion implosions due to shock-driven kinetic fuel-species stratification and viscous heating.  Physics of Plasmas, 25:056310, 2018.
  * [71] K. Molvig, N. M. Hoffman, B. J. Albright, E. M. Nelson, and R. B. Webster.  Knudsen Layer Reduction of Fusion Reactivity.  Physical Review Letters, 109:095001, 2012.
  * [72] B. J. Albright, K. Molvig, C.-K. Huang, A. N. Simakov, E. S. Dodd, N. M. Hoffman, G. Kagan, and P. F. Schmit.  Revised Knudsen-layer reduction of fusion reactivity.  Physics of Plasmas, 2013.
  * [73] G. Kagan, D. Svyatskiy, H. G. Rinderknecht, M. J. Rosenberg, A. B. Zylstra, C.-K. Huang, and C. J. McDevitt.  Self-similar structure and experimental signatures of suprathermal ion distribution in inertial confinement fusion implosions.  Physical Review Letters, 115:105002, 2015.