# Active turbulence and spontaneous phase separation in inhomogeneous
extensile active gels

Renato Assante  SUPA, School of Physics and Astronomy, The University of
Edinburgh, James Clerk Maxwell Building, Peter Guthrie Tait Road, Edinburgh,
EH9 3FD, United Kingdom  Dom Corbett  SUPA, School of Physics and Astronomy,
The University of Edinburgh, James Clerk Maxwell Building, Peter Guthrie Tait
Road, Edinburgh, EH9 3FD, United Kingdom  Davide Marenduzzo  SUPA, School of
Physics and Astronomy, The University of Edinburgh, James Clerk Maxwell
Building, Peter Guthrie Tait Road, Edinburgh, EH9 3FD, United Kingdom
Alexander Morozov  SUPA, School of Physics and Astronomy, The University of
Edinburgh, James Clerk Maxwell Building, Peter Guthrie Tait Road, Edinburgh,
EH9 3FD, United Kingdom

###### Abstract

We report numerical results for the hydrodynamics of inhomogeneous lyotropic
and extensile active nematic gels. By simulating the coupled Cahn–Hilliard,
Navier–Stokes and Beris–Edwards equation for the evolution of the density,
flow and orientational order of an active nematic, we ask whether density
variations are important to determine its emergent physics. As in constant-
density active gels, we find that increasing either activity or nematic
tendency (e.g., overall active matter density) triggers a transition between
an isotropic passive phase and an active nematic one. We show that density
inhomogeneities are important in the latter phase, where we find three types
of possible dynamical regimes. First, we observe regular patterns with defects
and vortices: these exist close to the passive-active transition. Second, for
larger activity, or deeper in the nematic phase, we find active turbulence, as
in constant-density active gels, but with exceedingly large density variation.
In the third regime, which is unique to inhomogeneous active nematics and
occurs for large nematic tendency and intermediate activity, we observe
spontaneous microphase separation into active and passive domains. The
microphase separated regime is notable in view of the absence of an explicit
demixing term in the underlying free energy which we use, and we provide a
theoretical analysis based on the common tangent construction which explains
its existence. We hope this regime can be probed experimentally in the future.

##  I Introduction

Active gels [1, 2, 3] are fascinating examples of non-equilibrium soft matter.
Some well-known realisations of these systems include solutions of
cytoskeletal filaments, such as actin or microtubules, interacting with
molecular motors, such as myosin or kinesin [4, 5]. Other instances come from
living materials, and encompass microbial suspensions of algae or bacteria [6,
7]. In an active gel, the constituent particles exert non-thermal forces on
their environments. Such forces can be modelled, at the simplest level, as
force dipoles, whose direction defines a nematic order parameter, which is a
fundamental quantity to describe the emergent physics of these systems [8].

The activity arising from the distribution of force dipoles leads to a
phenomenology which is strikingly different from that of passive colloidal
particles or polymer suspensions. For instance, these materials harbour a
“spontaneous flow” instability, which sets in for sufficiently strong
activity, and comprises a non-equilibrium transition between a quiescent
suspension and a state where activity fuels continuous motion [9, 10, 11, 12,
13]. For sufficiently large activity, the flow and orientation patterns of the
spontaneously flowing state are seemingly chaotic, and the associated state is
known as “active turbulence” [14, 15, 16, 17, 18, 19]. In active turbulence,
active gels self-organise into a random arrangement of vortices. Experiments
and theories suggest that in the nematic phase these vortices have a typical
length scale, arising from the competition between activity and elasticity
[14], while recent work points to important fundamental differences between
active turbulence and its more widely studied passive counterpart [18, 20].
Active gels also possess strongly non-Newtonian rheological properties [21],
such as marked activity-induced thinning or thickening [6, 11, 13, 22, 23],
Darcy-like flow [24], or negative drag in microrheology [22, 25].

Existing theories and simulations of active gel hydrodynamics typically
consider systems with constant density. In contrast, inspection of active
turbulent patterns found with microtubule-kinesin mixtures in the presence of
polyethylene-glycol (which causes adsorption to the oil-water interface [5,
26]) shows that the density of active material is significantly inhomogeneous.
While a linear stability analysis shows that in extensile gels, such as a
microtubule-kinesin mixture, the onset of spontaneous flow depends on
orientational bend fluctuations and density fluctuations should be irrelevant
[27], active turbulence is a highly non-linear phenomenon and the relevance of
density inhomogeneities to its physics remains unclear. Additionally, passive
colloidal particles aggregate in active nematics [28], through a mechanism
reminiscent of path coalescence [29, 30] or fluctuation-dominated phase
ordering [31]. This nonequilibrium aggregation shows that even a one-way
coupling between density and spontaneous flow (as tracers are affected by the
spontaneous flow but do not modify it) can in principle give rise to density
inhomogeneities, and more in general to nontrivial physics outside the reach
of a constant-density approximation.

To understand the role of density inhomogeneities in the physics of
spontaneous flow and active turbulence, here we study the hydrodynamic
equations of motion of a mixture of an isotropic fluid and an active nematic
gel by means of computer simulations. While the overall system is always
incompressible, the active gel component can change density, as is realistic
for the microtubule-kinesin mixtures considered in [26]. With respect to
previous work on Cahn–Hilliard models coupled to active nematics focussing on
the crossover between wet and dry systems driven by friction with the
substrate [32, 33, 34], here our focus is specifically on the qualitative role
of density inhomogeneities on the emerging physics and patterns. As in
constant-density active gels, we find that the system displays a transition
between a passive isotropic phase and an active nematic phase. This transition
can be triggered either by increasing activity or the nematic tendency of the
system (more specifically the coupling between active matter density and
orientational order). In the active nematic phase, though, density
inhomogenities play a fundamental role and give rise to some unique dynamical
behaviour. We find three dynamical regimes in this phase. Close to the
transition boundary, our lyotropic system settles into regular flowing
patterns, with approximately ordered spiral defect arrangements creating a
rotational active flow consisting of long-lived and stable vortices. For
larger activity, or deeper in the nematic phase, the flow becomes chaotic and
we observe active turbulence. Both regular patterns and active turbulence are
regimes which can be found in constant-density active gels. There are
differences though, as in our case the thermodynamic coupling between
orientational and compositional order parameters favours a density minimum, or
relative void of active matter, at defect cores. Additionally, active
turbulent patterns are characterised by a very broad distribution of local
densities. The last regime we observe is unique to inhomogeneous active
nematics, and is found even deeper in the nematic phase with respect to active
turbulence, but for low or intermediate activities. Here the active mixture
spontaneously phase separates into low-density disordered and high-density
nematic domains of irregular shape. We show by a semi-analytical theoretical
analysis that this phase separation is due to the coupling between density and
nematic ordering and hence is driven thermodynamically. The corresponding
coarsening is arrested in our case by the spontaneous active flow, and the
size of the steady state domains decreases with activity. We conclude by
discussing ways in which our study can be taken forward, both theoretically
and experimentally.

##  II Equations of motion

To describe the equilibrium properties of an inhomogeneous active nematic
system in the passive phase (i.e., when the activity parameter defined below
is switched off), we employ a Landau-de Gennes free energy ${\cal F}$, whose
density we call $f$. The latter consists of a contribution which characterises
the orientational order of the active liquid crystal, measured by a nematic
tensor $Q_{\alpha\beta}$, and of another contribution determining the physics
of density fluctuations, depending on the compositional order parameter $\phi$
(which measures the local density of active material). The liquid crystalline
free energy density we use is a standard one to describe passive nematic
liquid crystals [35], and is explicitly given by

$\displaystyle f_{\rm
LC}=\frac{A_{0}}{2}\left(1-\frac{\gamma(\phi)}{3}\right)Q_{\alpha\beta}^{2}-\frac{A_{0}\gamma(\phi)}{3}Q_{\alpha\beta}Q_{\beta\gamma}Q_{\gamma\alpha}$
$\displaystyle+\frac{A_{0}\gamma(\phi)}{4}\left(Q_{\alpha\beta}^{2}\right)^{2}+\frac{K}{2}\left(\partial_{\gamma}Q_{\alpha\beta}\right)^{2},$
(1)

where the first term is a bulk contribution, describing the isotropic-nematic
transition, while the second term is an elastic distortion term. In the
equation above, $A_{0}$ is a constant, $\gamma(\phi)$ controls the magnitude
of order (so that it may be viewed as an effective temperature or
concentration for thermotropic and lyotropic liquid crystals respectively),
while $K$ is an elastic constant — note we are using the (standard in this
field) one-constant approximation [35]. Here and in what follows Greek indices
denote cartesian components and summation over repeated indices is implied.
The coupling between concentration and ordering arises through $\gamma(\phi)$,
which equals

$\displaystyle\gamma(\phi)=\gamma_{0}+\phi({\bf r},t)\Delta,$ (2)

where $\gamma_{0}$ and $\Delta$ are appropriate constants. In our simulations
described below, we fix $\gamma_{0}$ and vary $\Delta$ (see Section 3 for full
parameter list).

The free energy density used to describe density fluctuations is instead given
by a simple function, used to describe binary fluid in the mixed (non-phase-
separating) regime, and its form is simply given by

$\displaystyle f_{\phi}=\frac{a}{2}\phi^{2},$ (3)

where $a$ is a constant related to the compressibility of the active gel
component. Note we do not include a surface-tension-like square gradient term,
$\frac{k_{\phi}}{2}\left(\partial_{\alpha}\phi\right)^{2}$, required for
stabilisation in conventional binary mixture models, as this is not necessary
in our study. Note also that the total free energy density $f=f_{\rm
LC}+f_{\phi}$ has two contributions which depend on $\phi$: besides the mixing
free energy $f_{\phi}$, Eq. 3, there is also a $\phi$ dependence in the liquid
crystalline free energy $f_{\rm LC}$, Eq. 1, through the $\gamma(\phi)$ term.

The fluid velocity, ${\bf u}$, obeys the continuity equation and the
Navier–Stokes equation,

$\displaystyle\qquad\qquad\partial_{\alpha}u_{\alpha}=0,$ (4)
$\displaystyle\rho\left(\partial_{t}+u_{\beta}\partial_{\beta}\right)u_{\alpha}=-\partial_{\alpha}p_{0}+\eta\partial_{\beta}^{2}u_{\alpha}$
$\displaystyle\qquad\qquad\qquad+\partial_{\beta}\Pi_{\alpha\beta}-\zeta\partial_{\beta}\left(\phi
Q_{\alpha\beta}\right),$ (5)

where $\rho$ is the fluid density, $\eta$ is an isotropic viscosity, and

$\displaystyle\Pi_{\alpha\beta}=2\xi(Q_{\alpha\beta}+{1\over
3}\delta_{\alpha\beta})Q_{\gamma\epsilon}H_{\gamma\epsilon}$ (6)
$\displaystyle-\xi H_{\alpha\gamma}(Q_{\gamma\beta}+{1\over
3}\delta_{\gamma\beta})-\xi(Q_{\alpha\gamma}+{1\over
3}\delta_{\alpha\gamma})H_{\gamma\beta}$
$\displaystyle-\partial_{\alpha}Q_{\gamma\nu}{\partial
f\over\partial\partial_{\beta}Q_{\gamma\nu}}+Q_{\alpha\gamma}H_{\gamma\beta}-H_{\alpha\gamma}Q_{\gamma\beta}.$
(7)

is the stress tensor, where $\zeta$ is the activity [8], and measures the
strength of active force dipoles. With the sign convention chosen here
$\zeta>0$ means extensile rods and $\zeta<0$ means contractile ones [8]. The
molecular field ${\bf H}$ which provides the driving motion is given by

${\bf H}=-{\delta{\cal F}\over\delta{\bf Q}}+({\bf I}/3){\rm Tr}{\delta{\cal
F}\over\delta{\bf Q}},$ (8)

where Tr denotes the tensorial trace.

The equation of motion for Q is taken to be [36]

$\displaystyle\partial_{t}{\bf Q}+{\bf u}\cdot{\nabla\bf Q}=\Gamma{\bf H}+{\bf
S},$ (9)

where $\Gamma$ is a collective rotational diffusion constant. The first term
on the left-hand side of Eq. (9) is the material derivative describing the
usual time dependence of a quantity advected by a fluid with velocity ${\bf
u}$. This is generalized for rod-like molecules by a second term

$\displaystyle{\bf S}=\left(\xi{\bf D}+{\bm{\omega}}\right)\cdot\left({\bf
Q}+\frac{1}{3}{\bf I}\right)+\left({\bf Q}+\frac{1}{3}{\bf
I}\right)\cdot\left(\xi{\bf D}-{\bm{\omega}}\right)$
$\displaystyle\qquad\qquad-2\xi\left({\bf Q}+\frac{1}{3}{\bf I}\right)\,{\rm
Tr}\left({\bf Q}\cdot\nabla{\bf u}\right),$ (10)

where

$\displaystyle{\bf D}=\frac{\nabla{\bf u}+\nabla{\bf u}^{\dagger}}{2},$ (11)
$\displaystyle{\bm{\omega}}=\frac{\nabla{\bf u}-\nabla{\bf u}^{\dagger}}{2}$

are the symmetric part and the anti-symmetric part respectively of the
velocity gradient tensor $\partial_{\beta}u_{\alpha}$. The constant $\xi$
depends on the molecular details of a given liquid crystal. The first term on
the right-hand side of Eq. (9) describes the relaxation of the order parameter
towards the minimum of the free energy.

Finally, the active material density, $\phi$, obeys a Cahn–Hilliard-like
equation,

$\displaystyle\partial_{t}{\bf\phi}+{\bf u}\cdot{\nabla\phi}=M\nabla^{2}\mu,$
(12)

where $\mu=\frac{\delta{\cal F}}{\delta\phi}$ is the chemical potential of the
active mixture, and $M$ is a mobility, which for simplicity we consider to be
constant.

##  III Numerical method

To study the dynamics of Eqs.(4)-(12), here we perform direct numerical
simulations. Note that in our simulations we assume that the fields are two-
dimensional and that the ${\bf Q}$ tensor describes nematic order in a 2D
plane (i.e., we assume that there is no out-of-plane nematic order). We employ
an in-house MPI-parallel code developed within the Dedalus spectral framework
[37]. Simulations are performed on a periodic rectangular domain
$[0,H]\times[0,H]$ with $H=200$ (here and below, all quantities are given in
simulation units). All fields are represented by de-aliased, double periodic
Fourier series with $512\times 512$ Fourier modes. Our time-iteration scheme
employs a 4th-order semi-implicit backward differentiation formula (BDF)
scheme [38] with the timestep $dt=0.1$. We have confirmed that our spatial and
temporal accuracy is sufficient to obtain numerically converged results. To
obtain statistically converged averages, simulations are performed for $30000$
time units.

In what follows, we fix $\rho=2$, $\eta=5/3$, $\xi=0.7$, $A_{0}=0.1$,
$K=0.01$, $\gamma_{0}=2$, $\Gamma=1$, and $a=0.003$, and $M=4$. These
parameters are chosen as they are in line with those used in previous hybrid
lattice Boltzmann simulations of both constant-density active nematics [11,
12] and self-motile active gel droplets [39].

##  IV Results

To discuss our results, we first present numerical simulations, and then a
semi-analytical discussion of phase separation in the passive limit which
sheds light on the phase diagram which we find.

###  IV.1 Simulation results

To address the role of density inhomogeneities in the hydrodynamics of
extensile active gels, we first report numerical simulations (for methodology
details, see Section 3). We characterise the dynamical behaviour of the system
as a function of two key parameters: activity, quantified by $\zeta$, and
tendency to acquire nematic order, quantified by
$\gamma(\phi_{0})=\gamma_{0}+\phi_{0}\Delta$. In practice, we change $\zeta$
and $\Delta$ in our simulations, keeping other parameters fixed (see Section 3
for a full list).

For each set of parameter values, we compute: (i) the average largest
eigenvalue of ${\bf Q}$, denoted by $\langle q\rangle$, which we use to
quantify the global magnitude of order; (ii) the average fluid velocity,
$\langle|{\bf u}|\rangle$; and (iii) the average variance of the compositional
order parameter. In each case, these averages are computed first spatially,
over a configuration, and then over time. We acquire data when the system is
in a statistical steady state, removing any initial transient. These
quantities allow us to build a phase diagram.

The magnitude of order, $\langle q\rangle$, is shown in Fig. 1 and allows us
to map the boundary between the isotropic and nematic phase. The nematic phase
can be reached by increasing either $\Delta$ (thereby $\gamma(\phi_{0})$) or
$\zeta$. The former is a thermodynamic route, the latter is a non-equilibrium
one. The non-equilibrium transition between isotropic and nematic phases
arises because activity-induced flows create shear, which in turn generates
nematic order. Analogous transitions in systems with _constant_ density have
already been studied [32, 40]. A linear stability calculation analogous to
that of [40] shows that an isotropic state of constant $\phi$ is unstable, for
sufficiently large $\zeta$, to an active nematic state with non-zero flow and
order. The critical threshold in the $\left(\Delta,\zeta\right)$ plane is
described by $\zeta_{\textnormal{c}}(\Delta)$, where

$\zeta_{\textnormal{c}}(\Delta)=\frac{2}{3}\xi
A_{0}\left(1-\frac{\gamma_{0}+\phi_{0}\Delta}{3}\right)\left(1+\frac{9\Gamma\eta}{2\xi^{2}}\right),$
(13)

corresponding to a straight line on the $(\Delta,\zeta)$ plane. This predicted
boundary is shown as a yellow line in Fig. 1. It agrees well with our numerics
for small $\Delta$, sufficiently far from the passive isotropic–nematic
transition (at $\zeta=0$ and $\Delta=0.7$). As the latter is a first-order
discontinuous transition, it is unsurprising that there is a quantitative
discrepancy with our linearised calculation close to this point.

Figure 1: Heatmap of the magnitude of nematic order as a function of $\zeta$
and $\Delta$. The plot can be used to define regions in parameter space where
the system is in the isotropic passive (bottom left) or active nematic (top
right) phase. The prediction for the onset of spontaneous flow (passive-active
transition) from linear stability analysis (see text) is shown as a yellow
line. Note that the isotropic–nematic transition in the passive limit
($\zeta=0$) occurs at $\gamma(\phi_{0})=2.7$ (corresponding to $\Delta=0.7$;
cf. section IV.2).

Throughout the nematic phase, except at $\zeta=0$, we find a non-zero flow in
steady state. In other words, the emerging nematic structures are always
active in the parameter range which we have explored. To characterise their
behaviour, we show in Fig. 2 the dynamical regimes we observe, before
discussing the phase diagram quantitatively.

Figure 2: From top to bottom rows represent: snapshots of density field,
$\left\\{\phi(x,y)\mid x,y\in[0,200]\right\\}$; director field, $\bm{\hat{n}}$
(scaled by local magnitude of order, $q$), overlaid on density (for
$x\in[87,113]$, $y\in[168,194]$); vorticity ($\omega_{z}$); largest eigenvalue
of $\bm{Q}$ tensor ($q$); $|Q_{xy}|$, qualitatively corresponding to a
Schlieren pattern as could be obtained experimentally with crossed polarisers.
Parameters are as in Section except: (a-e) $\zeta=0.16$, $\Delta=0.3$; (f-j)
$\zeta=0.2$, $\Delta=1$; (k-o) $\zeta=0.02$, $\Delta=1$. Colour scales are
shared between (f)-(j) and (k)-(o).

Fig. 2(a–e) shows an example of regular active patterns. These are found close
to the isotropic–nematic transition for moderate values of $\Delta$. The
example shown features regular density modulations in steady state (Fig. 2a),
where local minima are placed at spiral defects of topological charge $1$. The
latter are most easily visible through the 4-brush pattern in the Schlieren-
like plot of $|Q_{xy}|$ in Fig. 2e [35], and also correspond to deep minima in
the local nematic order parameter, as shown in Fig. 2d. There are also steady
vortices associated with the pattern, because spirals continuously rotate as
in previous models of defects in constant-density active nematics [41]. Such
vortices can be identified as maxima and minima in the vorticity plot in Fig.
2c. It is the thermodynamic coupling (proportional to $\Delta$) between order
parameter and density in our free energy that drives the local density
depletion at the centres of spirals, giving rise to a correlation between
density and order, where larger density is associated with larger order, and
vice versa. Our simulations show that the regular defect patterns we find
close to the transition can either be stationary or self-motile (see Suppl.
Movie 1 for the dynamics correspondent to the snapshot in Fig. 2a–e).

Fig. 2(f–j) shows an example of a different dynamical regime, obtained deeper
in the nematic phase, for sufficiently large $\zeta$ and $\Delta$ (see also
Suppl. Movie 2). Here, the patterns are never stable but display a chaotic
dynamics and diffuse around in the system. In line with constant-density
active nematic literature [15, 18, 20], we refer to these spatiotemporally
varying patterns as active turbulence. Our simulations show that active
turbulent patterns are accompanied by the appearance of defects — mostly of
half-integer topological charge, corresponding to two-brush patterns in
$|Q_{xy}|$ (Fig. 2j) — as in constant-density active gels. We also find this
regime is characterised by large density variations (Fig. 2f). As in the case
of regular patterns, here too the density field tends to decrease close to
defect cores (Fig. 2f,g). As the active chaotic spontaneous flow moves defects
around, streaks of density voids form (Fig. 2f,g) which follow defect
trajectories.

While the regimes in Figs. 2(a–e) and 2(f–j) are qualitatively reminiscent of
those found in the literature for constant-density active extensile gels [11,
32, 40], we also find a third dynamical regime which is unique to
inhomogeneous active nematics. This regime is found for low or intermediate
activity and sufficiently large $\Delta$, and consists in spontaneous phase
separation into active and passive domains (Fig. 2(k–o), see also Suppl. Movie
3). Active domains are nematically ordered, and their size depends on $\zeta$.
Importantly, this is an example of microphase separation, or arrested phase
separation, as coarsening does not proceed indefinitely and there are multiple
domains in steady state (Fig. 2k,l). In other words, the late-time domain size
does not scale with system size, and is instead determined by other
parameters: for instance, it decreases with increasing $\zeta$. The presence
of this spontaneous microphase separation regime is at first sight surprising,
as the $\phi$-dependent part of the free energy does not promote demixing by
itself. As we show in section IV.2, however, phase separation is driven by the
coupling between nematic order and local density of active matter.

To quantitatively delineate the phase diagram of the system, and the
boundaries between the three different dynamical regimes shown in Fig. 2, we
proceed as follows. We use the plot in Fig. 1 to identify the phases as
isotropic ($\langle q\rangle\simeq 0$) or nematic ($\langle q\rangle\neq 0$).
We then classify cases where the kinetic energy reaches a plateau as regular
patterns. To discriminate between active turbulence and spontaneous phase
separation, we look at probability distribution functions for $\phi$
(calculated over space, and averaged over time, see Fig. 3 for examples). If
this distribution is bimodal (i.e., it has two maxima), then we classify the
pattern as phase separated (Fig. 3, yellow curve). A similar identification
would arise from analysing the average density variance and the average flow
magnitude (Fig. 4). Flow and density variance are highest in the active
pattern and spontaneous phase separation regimes respectively. The phase
diagram for our inhomogeneous extensile mixture is shown in Fig. 5.

Figure 3: Probability distribution functions showing frequency versus local
density at a grid point, for the active pattern (blue curve), the active
turbulence (maroon curve) and the spontaneous phase separation (yellow curve)
regimes. It can be seen that the distribution corresponding to the active
turbulent regime is unimodal, while that associated with the spontaneous phase
separated one is bimodal, with peaks at $\phi\approx 0.19,1.41$, cf. binodal
points in Fig. 6. Figure 4: Combined heatmap of the averaged flow magnitude,
$|{\mathbf{u}}|$ (blue), and variance of $\phi$ (yellow), as functions of
$\zeta$ and $\Delta$. The plot can be used to identify the regions of
parameter space where phase separation is most prominent (upper left) or where
the steady active flow is strongest (right).

The density distribution plots show that density variations are a generic
feature throughout the active nematic phase. First we note that, without such
density variations, the spontaneous microphase separation regime could not
arise. The existence of this regime in a realisable system depends on the
value of $\Delta$, which will be determined by microscopic parameters.
Assuming an Onsager-like theory for density-dependent orientational order, we
expect this could be achievable with rod-like active particles of sufficiently
large aspect ratio. Note that we actually observe a bimodal distribution even
in the regular pattern regime, however here the magnitude of the peaks differs
by more than one order of magnitude. This corresponds to a detectable
depletion of active material at the cores of long-lived defects rather than
actual phase separation. A second observation is that the width of the density
distributions is remarkably large in the active turbulent regime, as the
sampled values of $\phi$ can vary from $\sim 0$ to $\sim 2$. This range is
similar to the spread observed in the spontaneous phase separation regime,
although the distribution remains peaked at $\phi\simeq\phi_{0}=1$ for active
turbulent patterns. This observation is in line with experimental microscopy
of quasi-2D extensile mixtures of microtubules and molecular motors, which
shows a highly inhomogeneous density of active material [26].

Figure 5: Phase diagram of the system in the $(\zeta,\Delta)$ plane. Each
square corresponds to one simulation. Regions A–D correspond respectively to
isotropic passive, active patterns, active turbulence, and spontaneous
microphase separation regimes.

###  IV.2 Coupling-induced demixing in the passive limit

To shed more light on the physics underlying the occurrence of our spontaneous
microphase separated regime, it is instructive to consider the passive limit
of our nematic mixture, where the behaviour is solely determined by
thermodynamic considerations. In this case, the free energy is minimised in
equilibrium, and a full phase diagram can be computed semi-analytically using
the common tangent construction (Fig. 6).

Figure 6: Coupling-induced demixing derived semi-analytically from the free
energy for a uniformly aligned passive nematic phase given in Eq. 14. Panel
(a) shows the free energy, $f_{\text{hom}}\left({q}_{0}(\phi),\phi\right)$,
which, for a given value of $\phi$, is minimized by ${q}_{0}(\phi)$ (see
Section IV.2). The second derivative with respect to $\phi$ is shown with a
dashed curve, and is negative in the spinodal region, highlighted in yellow.
Here the system is linearly unstable to demixing into regions of high and low
density, the values of which are given by the respective binodal points, shown
in grey. These points are found by the double tangent construction illustrated
by the green line, and share the same chemical potential and pressure ($\Pi$,
dashed in blue). In this plot $\Delta$ is fixed at 0.7, and the presence of
the lower spinodal point at $\phi=1$ is consistent with our simulation
results, as illustrated in Fig. 5. Panel (b) shows the $\phi$–$\Delta$ plane,
which is a phase diagram derived from the loci of the spinodal and binodal
points illustrated in panel (a), considered as we vary $\phi$ and $\Delta$.
The shaded spinodal and binodal regions are bounded by these loci, and show
the parameter combinations for which the system is unstable to phase
decomposition.

Recall that our model employs both conserved (composition $\phi$) and non-
conserved (nematic tensor $Q_{\alpha\beta}$) order parameters, which are
mutually coupled through the free energy given in Eq. 1. The resulting mixture
is an example of a lyotropic liquid crystal, such as that considered in [42],
albeit with an important distinction. Specifically, in the limit of no
coupling ($\Delta=0$ in Eq. 2), the system we consider exhibits no phase
separation and mixes freely. In what follows, therefore, we shall see that
passive phase separation is driven only by the coupling.

To find the equilibrium state, we consider a uniformly aligned (homogenous)
nematic phase, so that the elastic term contribution in Eq. 1 vanishes and the
total free energy can be written in terms of the magnitude of nematic order,
$q$, and of $\phi$ as follows,

$\displaystyle f_{\text{hom}}$
$\displaystyle=\frac{A_{0}}{3}\left(1-\frac{\gamma(\phi)}{3}\right)q^{2}-\frac{2A_{0}\gamma(\phi)}{27}q^{3}$
$\displaystyle\qquad\qquad+\frac{A_{0}\gamma(\phi)}{9}q^{4}+\frac{a}{2}\phi^{2}.$
(14)

We proceed from here with the assumption that our globally conserved
composition relaxes on a slower timescale than our non-conserved nematic
order, allowing us to find ${q}_{0}(\phi)$, the value of $q$ that minimizes
the total free energy for a given value of $\phi$, viz.

$\displaystyle{q}_{0}(\phi)$ $\displaystyle=\begin{cases}0,&\text{ if
}\gamma\leq\gamma_{c}\\\
\tfrac{1}{4}\left(1+\sqrt{9-24/\gamma(\phi)}\right),&\text{ if
}\gamma\geq\gamma_{c}\end{cases}$ (15)

We find that the resulting free energy,
$f_{\text{hom}}\left({q}_{0}(\phi),\phi\right)$, is a non-convex function of
$\phi$ for sufficiently large $\Delta$ (Fig. 6a). This pleasingly signals that
the nematic mixture is unstable to phase separation in the parameter range
(viz. ${\phi}_{0}=1$, $\Delta>0.7$) that matches our simulation results given
above.

In greater detail, the spinodal region — where the system is linearly unstable
to phase separation, for any perturbation however small — is precisely where
$f_{\text{hom}}\left({q}_{0}(\phi),\phi\right)$ is concave down. Thus it is
bounded by the inflexion points, where the second derivative of
$f_{\text{hom}}$ with respect to $\phi$,
$f_{\phi\phi}\left({q}_{0}(\phi),\phi\right)$, crosses the $\phi$ axis in Fig.
6a. The lower spinodal point is simply
${\phi}_{\text{c}}=(\gamma_{c}-\gamma_{0})/\Delta$, whereas the upper spinodal
point, ${\phi}_{\text{s}}$, was found numerically as the only real root of a
seventh order polynomial.

This effective free energy also gives, by a double tangent construction, the
binodal points, namely the values of the coexisting densities of the phase
separated systems. Specifically, as coexisting phases must share a common
chemical potential ($\mu=f_{\phi}\left({q}_{0}(\phi),\phi\right)$) and a
common pressure ($\Pi=f_{\text{hom}}\left({q}_{0}(\phi),\phi\right)-\mu\phi$),
we have a system of two equations in two unknowns. The two unknowns in
question are the two values of $\phi$ for which a straight line tangentially
touches the curve $f_{\text{hom}}\left({q}_{0}(\phi),\phi\right)$ exactly
twice, and were found numerically.

The loci of these spinodal and binodal points as we vary $\phi$ and $\Delta$
are the lines plotted in the $\phi$–$\Delta$ plane in Fig. 6b. Here we can see
that, in line with our simulations, there is a spinodal point at $\Delta=0.7$
for ${\phi}_{0}=1$ (compare with the phase diagram in Fig. 5).

We stress again that the analysis in this Section refers to the $\zeta=0$
passive limit of our mixture. In the active case, macroscopic phase separation
is generically arrested by the activity-induced spontaneous flow. This leads
to microphase separation at steady state, similarly to what has been observed
in sheared passive binary fluids [43] or in active model H [39].

Quantitatively, our theory predicts a spinodal point at $\Delta=0.7$ for
${\phi}_{0}=1$, which is consistent with our simulation results (Fig. 5).
Likewise, the binodal points for $\Delta=1$, viz. $\phi=0.18$ and $\phi=1.79$,
are consistent with the histogram peaks shown in Fig. 3, viz. $\phi=0.19$ and
$\phi=1.41$ — the lower value of the latter peak is likely due to the active
flow which drives the system away from thermodynamic equilibrium (assumed in
Fig. 6).

##  V Discussion and Conclusions

In summary, we have used computer simulations to study the hydrodynamics of an
inhomogeneous active nematic gel. With respect to conventional models for
active gels, which only consider the velocity field and ${\bf Q}$ tensor, our
theory also allows for the time evolution of the active matter density $\phi$.
Previous work has shown by a linear stability analysis that density
fluctuations are irrelevant for the physics of the “generic instability” of
active gels [27], which stands for the transition between the passive
(quiescent) and the active (spontaneously flowing) phase. It has however
remained unclear what their role is deep in the active phase, where
nonlinearities are important; shedding light on this issue has been the focus
of our current work.

Our main result is the quantitative characterisation of the phase diagram of
inhomogeneous active nematic (Fig. 5). We have found that there are three
regimes with distinct emergent behaviour in the active phase. First, close to
the transition between the passive isotropic and active nematic phase, there
are regular patterns typically composed of self-assembled rotating spirals.
Second, deeper in the active nematic phase there is an active turbulent regime
featuring chaotic dynamics of vortices and half-integer nematic defects.
Third, for low or moderate activity and large enough nematic tendency
($\Delta$ in our phase diagram in Fig. 5), we find spontaneous phase
separation into active and passive domains. This latter phase separation is
arrested by the active flow, so that domains do not coarsen past a typical
size, which decreases with increasing activity.

The regular spiral/vortex patterns we find are reminiscent of those observed
with polar active gels in the ordered phase [41, 44]. While polar nematics can
only exhibit defects with integer topological charge, defects in active
(apolar) nematics normally have half-integer topological charge, so it is non-
trivial that close to the passive-active transitions we observe spirals, whose
topological charge is $+1$. 111Note that in our geometry integration of the
topological charge density, defined as in [52], is conserved and equal to $0$.
Notwithstanding this qualitative resemblance, the pattern we observe also have
a non-trivial spatiotemporal dynamics (Suppl. Movie 1). The core of our
spirals are also associated with notable density minima, or voids, which arise
because of the coupling between nematic order and density in the free energy
of the system. The mechanism responsible for this coupling is the same that
drives inert colloidal particles or isotropic droplets (with no anchoring on
their surface) to the defect cores or disclinations in passive liquid crystals
[46, 28].

The chaotic, active turbulent regime we find for sufficiently large $\zeta$ is
an analog of the homonymous regime in constant-density active gels [47, 48,
18, 20]. An important feature of this regime in our simulations, though, is
that there are very large density fluctuations (Fig. 3a). These are
qualitatively in line with experimental observations of active turbulence in
microtubule-motor mixtures, which show substantial inhomogeneities in
microtubule density over a sample [5, 26]. Whether such density variations
lead to a fundamental change in the scaling properties of active turbulence is
an open question which we believe deserves further investigation, for instance
by a quantitative detailed analysis of the scaling of velocity-velocity
correlations [16, 17, 49].

Regarding the spontaneous microphase separated regime, this is notable
especially because there is no term in the free energy density $f_{\phi}$
which explicitly favours demixing. In other words, in the absence of $f_{\rm
LC}$ the system would remain uniform. Phase separation arises due to the
coupling between density and order, measured by the parameter $\Delta$. In
this sense, phase separation is not driven by activity but rather
thermodynamically, and indeed it can be explained with a theoretical
discussion of the free energy in the passive limit (Fig. 6). In simulations we
observe a microphase separated pattern rather than macroscopic phase
separation, as the active flow arrests coarsening, and controls the size of
the steady-state domains observed at late times, similarly to the case of
active model H [39]. While experimental realisation of active nematics have
shown plenty of instances of active turbulence [5, 26], the spontaneous
microphase separation regime appears to not have been found in the lab yet.
Our model suggests that the most promising avenue to realise this regime
experimentally is to control the density-order coupling $\Delta$. The latter
may be estimated by monitoring how the isotropic-nematic transition point
depends on the density of nematogenic particles (for instance, microtubules)
in the passive limit of no activity.

We close our result discussion with a note of caution. While the density of
active matter varies, the overall fluid is incompressible: in other words, the
total density of active matter and underlying solvent is constant.
Incompressibility is important especially for high activity and intermediate
$\Delta$ — i.e., in the bottom right region in the phase diagram in Fig. 5,
which is challenging to characterise accurately. There we typically find
active turbulent patterns with vortex correlation length similar to the system
size. We found that these are replaced by an apparent phase separation into
high and low-density patches if the Navier–Stokes solver allows for fluid
compressibility and the Mach number is not sufficiently small, which occurs
when using, for instance, hybrid Lattice Boltzmann simulations at $\zeta>0.1$.
These states are artifacts in the current model, and using a purely
incompressible and 2D fluid they are not found, in line with simulations with
tracers in turbulent flows which only show aggregation in compressible fluids,
or due to inertial effects [50, 31]. It would however be of interest to see
whether these patters may be recovered in thin active matter films, which can
behave effectively as compressible fluids [51].

Looking ahead, we can suggest a few directions in which our work can be
carried forward. First, it would be of interest to understand from a more
fundamental point of view the universal properties of the dynamical regimes we
have identified. For instance, one could quantify the dependence of vortex
correlation length and pattern size on physical parameters, and the scaling of
the power spectra of the kinetic energy. This would allow to clarify the
important theoretical question of whether or not inhomogeneous active
turbulence is in the same universality class of turbulence in constant-density
active gels. Second, from the experimental point of view, it would be
desirable to compare more quantitatively density distributions in active
turbulence with those predicted by our simulations. Third, with regards to
computer simulations, it would be highly interesting to explore the phase
behaviour and dynamics of inhomogeneous active nematics in 3D, comparing and
contrasting it with the one found here in 2D.

##  VI Acknowledgement

We thank S. Samatas for useful discussions. A. M. acknowledges support by the
Engineering and Physical Sciences Research Council (grant number
EP/V048198/1). For the purpose of open access, the author has applied a
Creative Commons Attribution (CC BY) licence to any Author Accepted Manuscript
version arising from this submission.

## References

  * Ramaswamy [2010] S. Ramaswamy, The mechanics and statistics of active matter, Annu. Rev. Condens. Matter Phys. 1, 323 (2010).
  * Prost  _et al._ [2015] J. Prost, F. Jülicher, and J.-F. Joanny, Active gel physics, Nat. Phys. 11, 111 (2015).
  * Doostmohammadi  _et al._ [2018] A. Doostmohammadi, J. Ignés-Mullol, J. M. Yeomans, and F. Sagués, Active nematics, Nat. Commun. 9, 1 (2018).
  * MacKintosh and Schmidt [2010] F. C. MacKintosh and C. F. Schmidt, Active cellular materials, Curr. Opin. Cell Biol. 22, 29 (2010).
  * Sanchez  _et al._ [2012] T. Sanchez, D. T. Chen, S. J. DeCamp, M. Heymann, and Z. Dogic, Spontaneous motion in hierarchically assembled active matter, Nature 491, 431 (2012).
  * Hatwalne  _et al._ [2004] Y. Hatwalne, S. Ramaswamy, M. Rao, and R. A. Simha, Rheology of active-particle suspensions, Phys. Rev. Lett. 92, 118101 (2004).
  * Wioland  _et al._ [2013] H. Wioland, F. G. Woodhouse, J. Dunkel, J. O. Kessler, and R. E. Goldstein, Confinement stabilizes a bacterial suspension into a spiral vortex, Phys. Rev. Lett. 110, 268102 (2013).
  * Marchetti  _et al._ [2013] M. C. Marchetti, J.-F. Joanny, S. Ramaswamy, T. B. Liverpool, J. Prost, M. Rao, and R. A. Simha, Hydrodynamics of soft active matter, Rev. Mod. Phys. 85, 1143 (2013).
  * Simha and Ramaswamy [2002] R. A. Simha and S. Ramaswamy, Hydrodynamic fluctuations and instabilities in ordered suspensions of self-propelled particles, Phys. Rev. Lett. 89, 058101 (2002).
  * Voituriez  _et al._ [2005] R. Voituriez, J.-F. Joanny, and J. Prost, Spontaneous flow transition in active polar gels, EPL 70, 404 (2005).
  * Marenduzzo  _et al._ [2007a] D. Marenduzzo, E. Orlandini, and J. Yeomans, Hydrodynamics and rheology of active liquid crystals: a numerical investigation, Phys. Rev. Lett. 98, 118102 (2007a).
  * Marenduzzo  _et al._ [2007b] D. Marenduzzo, E. Orlandini, M. Cates, and J. Yeomans, Steady-state hydrodynamic instabilities of active liquid crystals: Hybrid lattice boltzmann simulations, Phys. Rev. E 76, 031921 (2007b).
  * Cates  _et al._ [2008] M. Cates, S. Fielding, D. Marenduzzo, E. Orlandini, and J. Yeomans, Shearing active gels close to the isotropic-nematic transition, Phys. Rev. Lett. 101, 068102 (2008).
  * Giomi [2015] L. Giomi, Geometry and topology of turbulence in active nematics, Phys. Rev. X 5, 031003 (2015).
  * Doostmohammadi  _et al._ [2017] A. Doostmohammadi, T. N. Shendruk, K. Thijssen, and J. M. Yeomans, Onset of meso-scale turbulence in active nematics, Nat. Commun. 8, 1 (2017).
  * Stenhammar  _et al._ [2017] J. Stenhammar, C. Nardini, R. W. Nash, D. Marenduzzo, and A. Morozov, Role of correlations in the collective behavior of microswimmer suspensions, Phys. Rev. Lett. 119, 028005 (2017).
  * Bárdfalvy  _et al._ [2019] D. Bárdfalvy, H. Nordanger, C. Nardini, A. Morozov, and J. Stenhammar, Particle-resolved lattice boltzmann simulations of 3-dimensional active turbulence, Soft Matter 15, 7747 (2019).
  * Carenza  _et al._ [2020] L. N. Carenza, L. Biferale, and G. Gonnella, Cascade or not cascade? energy transfer and elastic effects in active nematics, EPL 132, 44003 (2020).
  * Kozhukhov and Shendruk [2022] T. Kozhukhov and T. N. Shendruk, Mesoscopic simulations of active nematics, Science Advances 8, eabo5788 (2022).
  * Alert  _et al._ [2020] R. Alert, J.-F. Joanny, and J. Casademunt, Universal scaling of active nematic turbulence, Nat. Phys. 16, 682 (2020).
  * Saintillan [2018] D. Saintillan, Rheology of active fluids, Annu. Rev. Fluid Mech. 50, 563 (2018).
  * Foffano  _et al._ [2012a] G. Foffano, J. Lintuvuori, A. Morozov, K. Stratford, M. Cates, and D. Marenduzzo, Bulk rheology and microrheology of active fluids, Eur. Phys. J. E 35, 1 (2012a).
  * Martinez  _et al._ [2020] V. A. Martinez, E. Clément, J. Arlt, C. Douarche, A. Dawson, J. Schwarz-Linek, A. K. Creppy, V. Škultéty, A. N. Morozov, H. Auradou, _et al._ , A combined rheometry and imaging study of viscosity reduction in bacterial suspensions, Proc. Natl. Acad. Sci. USA 117, 2326 (2020).
  * Mackay  _et al._ [2020] F. Mackay, J. Toner, A. Morozov, and D. Marenduzzo, Darcy’s law without friction in active nematic rheology, Phys. Rev. Lett. 124, 187801 (2020).
  * Foffano  _et al._ [2012b] G. Foffano, J. Lintuvuori, K. Stratford, M. Cates, and D. Marenduzzo, Colloids in active fluids: Anomalous microrheology and negative drag, Phys. Rev. Lett. 109, 028103 (2012b).
  * Martínez-Prat  _et al._ [2019] B. Martínez-Prat, J. Ignés-Mullol, J. Casademunt, and F. Sagués, Selection mechanism at the onset of active turbulence, Nat. Phys. 15, 362 (2019).
  * Baskaran and Marchetti [2009] A. Baskaran and M. C. Marchetti, Statistical mechanics and hydrodynamics of bacterial suspensions, Proc. Natl. Acad. Sci. USA 106, 15567 (2009).
  * Foffano  _et al._ [2019] G. Foffano, J. S. Lintuvuori, K. Stratford, M. Cates, and D. Marenduzzo, Dynamic clustering and re-dispersion in concentrated colloid-active gel composites, Soft Matter 15, 6896 (2019).
  * Maritan and Banavar [1994] A. Maritan and J. R. Banavar, Chaos, noise, and synchronization, Phys. Rev. Lett. 72, 1451 (1994).
  * Wilkinson and Mehlig [2003] M. Wilkinson and B. Mehlig, Path coalescence transition and its applications, Phys. Rev. E 68, 040101 (2003).
  * Das and Barma [2000] D. Das and M. Barma, Particles sliding on a fluctuating surface: phase separation and power laws, Phys. Rev. Lett. 85, 1602 (2000).
  * Thampi  _et al._ [2015] S. P. Thampi, A. Doostmohammadi, R. Golestanian, and J. M. Yeomans, Intrinsic free energy in active nematics, EPL 112, 28004 (2015).
  * Doostmohammadi  _et al._ [2016] A. Doostmohammadi, M. F. Adamer, S. P. Thampi, and J. M. Yeomans, Stabilization of active matter by flow-vortex lattices and defect ordering, Nat. Commun. 7, 1 (2016).
  * Thijssen  _et al._ [2021] K. Thijssen, D. A. Khaladj, S. A. Aghvami, M. A. Gharbi, S. Fraden, J. M. Yeomans, L. S. Hirst, and T. N. Shendruk, Submersed micropatterned structures control active nematic flow, topology, and concentration, Proceedings of the National Academy of Sciences 118, e2106038118 (2021).
  * De Gennes and Prost [1993] P.-G. De Gennes and J. Prost,  _The physics of liquid crystals_ , 83 (Oxford University Press, 1993).
  * Beris  _et al._ [1994] A. N. Beris, B. J. Edwards, _et al._ ,  _Thermodynamics of flowing systems: with internal microstructure_ , 36 (Oxford University Press, 1994).
  * Burns  _et al._ [2020] K. J. Burns, G. M. Vasil, J. S. Oishi, D. Lecoanet, and B. P. Brown, Dedalus: A flexible framework for numerical simulations with spectral methods, Phys. Rev. Res. 2, 023068 (2020).
  * Wang and Ruuth [2008] D. Wang and S. J. Ruuth, Variable step-size implicit-explicit linear multistep methods for time-dependent partial differential equations, J. Comput. Math. 26, 838 (2008).
  * Tiribocchi  _et al._ [2015] A. Tiribocchi, R. Wittkowski, D. Marenduzzo, and M. E. Cates, Active model h: scalar active matter in a momentum-conserving fluid, Phys. Rev. Lett. 115, 188302 (2015).
  * Santhosh  _et al._ [2020] S. Santhosh, M. R. Nejad, A. Doostmohammadi, J. M. Yeomans, and S. P. Thampi, Activity induced nematic order in isotropic liquid crystals, J. Stat. Phys. , 1 (2020).
  * Kruse  _et al._ [2004] K. Kruse, J.-F. Joanny, F. Jülicher, J. Prost, and K. Sekimoto, Asters, vortices, and rotating spirals in active gels of polar filaments, Phys. Rev. Lett. 92, 078101 (2004).
  * Matsuyama  _et al._ [2002] A. Matsuyama, R. Evans, and M. Cates, Non-uniformities in polymer/liquid crystal mixtures, Eur. Phys. J. E 9, 79 (2002).
  * Stansell  _et al._ [2006] P. Stansell, K. Stratford, J.-C. Desplat, R. Adhikari, and M. Cates, Nonequilibrium steady states in sheared binary fluids, Phys. Rev. Lett. 96, 085701 (2006).
  * Elgeti  _et al._ [2011] J. Elgeti, M. Cates, and D. Marenduzzo, Defect hydrodynamics in 2d polar active fluids, Soft Matter 7, 3177 (2011).
  * Note [1] Note that in our geometry integration of the topological charge density, defined as in [52], is conserved and equal to $0$.
  * Ravnik  _et al._ [2011] M. Ravnik, G. P. Alexander, J. M. Yeomans, and S. Žumer, Three-dimensional colloidal crystals in liquid crystalline blue phases, Proc. Natl. Acad. Sci. USA 108, 5188 (2011).
  * Thampi  _et al._ [2013] S. P. Thampi, R. Golestanian, and J. M. Yeomans, Velocity correlations in an active nematic, Phys. Rev. Lett. 111, 118101 (2013).
  * Linkmann  _et al._ [2019] M. Linkmann, G. Boffetta, M. C. Marchetti, and B. Eckhardt, Phase transition to large scale coherent structures in two-dimensional active matter turbulence, Phys. Rev. Lett. 122, 214503 (2019).
  * Škultéty  _et al._ [2020] V. Škultéty, C. Nardini, J. Stenhammar, D. Marenduzzo, and A. Morozov, Swimming suppresses correlations in dilute suspensions of pusher microorganisms, Phys. Rev. X 10, 031059 (2020).
  * Falkovich  _et al._ [2001] G. Falkovich, K. Gawȩdzki, and M. Vergassola, Particles and fields in fluid turbulence, Rev. Mod. Phys. 73, 913 (2001).
  * Voituriez  _et al._ [2006] R. Voituriez, J. Joanny, and J. Prost, Generic phase diagram of active polar films, Phys. Rev. Lett. 96, 028102 (2006).
  * Blow  _et al._ [2014] M. L. Blow, S. P. Thampi, and J. M. Yeomans, Biphasic, lyotropic, active nematics, Phys. Rev. Lett. 113, 248303 (2014).


# Mid-T Heat Treatments on BCP’ed Coaxial Cavities at TRIUMF

P. Kolb<EMAIL_ADDRESS>TRIUMF, 4004 Wesbrook Mall, Vancouver, B.C. V6T 2A3,
Canada  Z. Yao  TRIUMF, 4004 Wesbrook Mall, Vancouver, B.C. V6T 2A3, Canada
A. Blackburn  University of Victoria, Victoria, B.C., Canada  R. Gregory
University of Victoria, Victoria, B.C., Canada  TRIUMF, 4004 Wesbrook Mall,
Vancouver, B.C. V6T 2A3, Canada  D. Hedji  University of Victoria, Victoria,
B.C., Canada  TRIUMF, 4004 Wesbrook Mall, Vancouver, B.C. V6T 2A3, Canada  M.
McMullin  University of Victoria, Victoria, B.C., Canada  TRIUMF, 4004
Wesbrook Mall, Vancouver, B.C. V6T 2A3, Canada  T. Junginger  University of
Victoria, Victoria, B.C., Canada  TRIUMF, 4004 Wesbrook Mall, Vancouver, B.C.
V6T 2A3, Canada  R.E. Laxdal  TRIUMF, 4004 Wesbrook Mall, Vancouver, B.C. V6T
2A3, Canada

###### Abstract

Mid-T heat treatments in the range from 250 to 400 ∘C on superconducting
radio-frequency (SRF) cavities have been shown to provide high quality factors
that rise with applied rf field strength in high frequency, electro-polished
(EP), elliptical cavities operating at 2 K, similar to nitrogen doped
cavities. The rise in quality factor is attributed to a decrease in the
temperature dependent part of the surface resistance R${}_{\textrm{\scriptsize
BCS}}$. Until now, no results have been reported for these new treatments on
quarter-wave resonators (QWR) and half-wave resonators (HWR). The TRIUMF
multi-mode coaxial cavities are dedicated test cavities that allow frequency
and temperature resolved performance characterization of treatments without
changing environments, therefore providing an excellent test vehicle to test
these new treatments with rf frequencies ranging from 200 to 1200 MHz.

In this paper, performance measurements from both QWR and HWR cavities are
reported and their performance compared with four different treatments:
baseline, a conventional 120 ∘C low temperature bake for 48 hours, and two
mid-T bakes at 300 and 400 ∘C for 3 hours. In addition, sample analysis using
SEM, EDX and SIMS of witness samples is also shown.

It is found that the mid-T bakes are not directly transferable to low
frequency cavities. In the fundamental modes of the two test cavities, no
performance gain over the baseline treatment nor a decreasing temperature
dependent component with rising rf amplitude was observed. At frequencies
above 1 GHz and low temperatures, the mid-T bakes show a reduced field
dependence of R${}_{\textrm{\scriptsize BCS}}$ compared to both the baseline
and 120 ∘C treatments.

##  I Introduction

Modern nuclear physics experiments, such as the experiments conducted at the
ISAC-II [1] facility at TRIUMF, often rely on heavy ion particle accelerators
for high energy beams. To deliver these beams in an efficient way often
superconducting radio-frequency (SRF) technology is used. A major cost driver
for SRF based accelerators are cryogenic losses in SRF cavity walls caused by
the rf field used to accelerate the beam. The required cryogenic
infrastructure and liquid helium (LHe) as coolant are costly on capital and
operational budgets. Cryogenic losses are inversely proportional to the
unloaded quality factor Q${}_{\textrm{\scriptsize 0}}$. Therefore, increasing
the quality factor in SRF cavities is essential to reduce the cryogenic load
and therefore reduce operating costs or design accelerators with smaller
cryoplants. This is most relevant for SRF based accelerators with a high duty
factor, where the active load dominates over the static load. Additionally,
increasing the operational accelerating fields, and therefore usable voltage,
can lead to shorter, more compact accelerators, reducing real estate expenses.

Several novel heat treatments that increase Q${}_{\textrm{\scriptsize 0}}$
more than conventional treatments have been been developed and applied to 1.3
GHz single cell elliptical cavities [2, 3, 4], all on electro-polished (EP)
cavities. Coaxial style cavities such as quarter-wave and half-wave resonators
(QWR and HWR respectively) used in heavy ion accelerators are often treated
with buffered chemical polishing (BCP) due to the more complex geometry and
the resulting technical challenges of using EP.

Figure 1: Photo of the HWR (left) and QWR (right) cavities used in this study.

Recently discovered heat treatments at 250 to 500 ∘C in elliptical 1.3 GHz
cavities with an EP surface finish have been shown to produce comparable
performance to nitrogen doped cavities [5, 6, 7], the current state-of-the-art
treatment to date for high Q${}_{\textrm{\scriptsize 0}}$ at moderate
gradients [8], while being a simpler cavity processing recipe. The performance
increase is a result of oxygen dissolution and diffusion of the native niobium
oxide layers [9]. At temperatures between 250 and 500 ∘C the
Nb${}_{\textrm{\scriptsize 2}}$O${}_{\textrm{\scriptsize 5}}$ layer dissolves
and regrows back shallower after air exposure [10]. This avoids using an
external doping agent like nitrogen.

Figure 2: Temperature and residual gas profiles during the heat treatment of
the QWR at 300 ∘C. The other mid-T treatments follow similar patterns. The
fluctuations around goal temperature are due to the slow response of the
cavity temperature regulation.

These mid-T or ’oxygen diffusion’ treatments allow quicker processing compared
to N${}_{\textrm{\scriptsize 2}}$-doping as they eliminate the injection of
N${}_{\textrm{\scriptsize 2}}$ into the furnace at high temperatures and a
post-doping EP step. This simplified process decreases processing time over
N${}_{\textrm{\scriptsize 2}}$-doping and should also increase consistency in
larger cavity productions. Since mid-T bakes show a decreasing
R${}_{\textrm{\scriptsize BCS}}$ with increasing rf field amplitude, this
treatment is of high interest for low frequency SRF cavities that typically
operate at a temperature of 4.2 K. At this operating temperature, the surface
resistance R${}_{\textrm{\scriptsize s}}$ is dominated by its temperature
dependent component R${}_{\textrm{\scriptsize BCS}}$ and the temperature
independent part R${}_{\textrm{\scriptsize Res}}$ contributes only a small
fraction. A decreasing R${}_{\textrm{\scriptsize BCS}}$ at 4 K would be an
enormous benefit for SRF based heavy ion accelerators operating in continuous
wave (cw) mode at moderate peak surface fields, such as the ISAC-II or FRIB
accelerators. The first reported results on BCP’ed 1.3 GHz nine-cell cavities
[11] do not show the characteristic increase in Q${}_{\textrm{\scriptsize
0}}$, but still a larger Q${}_{\textrm{\scriptsize 0}}$ than a baseline
treatment. This is similar to the difference between EP’ed and BCP’ed
N${}_{\textrm{\scriptsize 2}}$-doped cavities [12]. Similarly, results on
EP’ed 650 MHz ellitpical cavities with either N${}_{\textrm{\scriptsize
2}}$-doping [13] or a 300 ∘C bake [14] do not show the increasing
Q${}_{\textrm{\scriptsize 0}}$, suggesting that there is a frequency
dependence to both treatments.

In this paper, the effects of two mid-T bakes at 300 and 400 ∘C on a QWR and a
HWR, both with BCP surface finishes, are studied and presented. The cavities
are shown in Fig. 1.

This paper presents the used methodology in Sec. II, results of cavity
measurements and sample studies in Sec. III and findings are discussed in Sec.
IV before summarizing them in Sec. V.

##  II Methodology

Two dedicated research cavities are used in this study: one QWR and one HWR,
both shown in Fig. 1. The two cavities and many aspects of the data
acquisition and processing are described in detail in [15]. These cavities are
simplified versions of the QWR and HWR geometries and designed for excitation
of not only the fundamental rf mode, but also higher order TEM modes. This
allows the extraction of the frequency dependence of key performance
parameters such as quench field and Q${}_{\textrm{\scriptsize 0}}$, defined as

$\displaystyle Q_{0}=\frac{\omega U}{P}$ (1)

with $\omega$ as resonant frequency, $U$ as stored energy in the RF field, and
$P$ as dissipated power in the cavity walls, without changing to different
cavities or different testing parameters such as cooldown speed or external
magnetic fields. These two cavities cover a frequency range from 220 MHz to
1.2 GHz, covering commonly used frequencies for SRF cavities.

Figure 3: (a) At 4 K the 220 MHz mode of the QWR with the 120 ∘C bake performs
significantly better with a lower surface resistance over a large range field
amplitudes up to the quench field. (b) At 2 K the mid-T bakes increase
G/Q${}_{\textrm{\scriptsize 0}}$ compared to the baseline and 120 ∘C bake.

This work compares performance of the two test cavities with a baseline
treatment (bulk BCP - 800 ∘C degassing - flash BCP), a commonly used low
temperature bake at 120 ∘C for 48 hours, and two mid temperature bakes at 300
and 400 ∘C for 3 hours. Between the 120 ∘C, 400 ∘C, and 300 ∘C heat treatments
the surface is reset with a 20 $\mu$m etch via BCP. Measurements of
Q${}_{\textrm{\scriptsize 0}}$ as a function of peak magnetic surface field
B${}_{\textrm{\scriptsize P}}$ are done at 4.2 K and 2.0 K, as well as
measurements of Q${}_{\textrm{\scriptsize 0}}$ during the cooldown between the
two temperatures at selected field amplitudes. This allows for the
deconvolution of the surface resistance R${}_{\textrm{\scriptsize s}}$ into
its temperature dependent and independent components,
R${}_{\textrm{\scriptsize BCS}}$ and R${}_{\textrm{\scriptsize Res}}$
respectively, by fitting the data to

$\displaystyle R_{s}(T)$
$\displaystyle=\frac{a_{0}}{T}\ln\left(\frac{4k_{b}T}{\hbar\omega}\right)\exp\left(-\frac{a_{1}(T)T_{c}}{T}\right)+a_{2}$
(2) $\displaystyle=R_{BCS}(T)+R_{Res}$ (3)

which is a parameterized version of the BCS surface resistance [16] with
$a_{0}$, $a_{1}$ and $a_{2}$ as fit parameters. The energy gap $\Delta$ is
represented by $a_{1}$ and its temperature dependence is modelled [17] as

$\displaystyle\frac{a_{1}(T)}{a_{1}(0\,K)}=\sqrt{\cos\left(\frac{\pi}{2}\left(\frac{T}{T_{c}}\right)^{2}\right)}.$
(4)

Further details of the test setup and data processing, including accounting
for the non-uniform field distribution across the surface [18], are described
in [15]. Measurements of Q${}_{\textrm{\scriptsize 0}}$ and
B${}_{\textrm{\scriptsize p}}$ are done in a zeroed external magnetic field
using a set of Helmholtz coils which surround the cavity and cryostat. In this
paper, the performance is often presented as G/Q${}_{\textrm{\scriptsize 0}}$
with G as the geometry factor defined as

$\displaystyle G=\frac{\omega\mu_{0}\int_{V}|H^{2}|dV}{\int_{S}|H^{2}|dS}$ (5)

while R${}_{\textrm{\scriptsize s}}$, R${}_{\textrm{\scriptsize Res}}$ and
R${}_{\textrm{\scriptsize BCS}}$ account for the non-uniform field
distribution. Not accounted for in the measurements are differences between
the helium bath temperature and the temperature on the rf surface, which can
lead to a measurable effect in Q${}_{\textrm{\scriptsize 0}}$ at higher field
amplitudes if the bath temperature is above the $\lambda$ point of LHe at 2.17
K [19].

The high and mid-T heat treatments are processed in TRIUMF’s induction
furnace, which is also described in [15] and [20]. The temperature is
controlled via an on/off controller set to the desired temperature. Due to the
temperature measurement and regulation setup, there is an oscillation around
the goal temperature as the measured temperature lags behind the temperature
of the heating element. A typical mid-T bake is shown in Fig. 2, showing the
heat treatment of the QWR at 300 ∘C for 3 hours. A residual gas analyser (RGA)
running during the heat treatment shows the partial pressures of selected
atomic masses, which correspond to molecules and elements such as
H${}_{\textrm{\scriptsize 2}}$, H${}_{\textrm{\scriptsize 2}}$O,
CO${}_{\textrm{\scriptsize 2}}$, N${}_{\textrm{\scriptsize 2}}$/CO and Ar. It
is important to note that the cavity ports are covered with caps made out of
Nb. The purpose of these caps is to prevent line of sight from potentially
contaminated furnace surfaces [21]. The caps are installed in a class 100
clean room environment after the cavity has been rinsed with high pressure
pure water as would be typically done before a clean assembly. The cavity is
then bagged and transported to the furnace and only just before insertion into
the heating chamber the bag is removed. During the heat treatment, a number of
witness samples were placed with the cavity in the furnace, but outside the rf
volume of the cavity. These samples are used for analysis with secondary
electron microscopy (SEM), energy-dispersive x-ray spectroscopy (EDX),
performed at the University of Victoria, and secondary ion mass spectroscopy
(SIMS), measured at the University of Western Ontario. All cavity cryogenic
measurements were conducted at TRIUMF’s SRF facilities.

Figure 4:  (a) 4 K measurement of the 390 MHz mode of the HWR show similar
behaviour to the QWR fundamental mode, fig 3, with the 120 ∘C bake out
performing all others. (b) The baseline and 120 ∘C bake perform similarly,
while the mid-T bakes add to R${}_{\textrm{\scriptsize s}}$. Figure 5:  (a) At
644 MHz and 4 K first differences between the mid-T bakes and the baseline
become noticeable, (b) but at 2 K both the baseline and 120 ∘C still produce
the highest Q${}_{\textrm{\scriptsize 0}}$ at low field.

The 120 ∘C bake [22] has previously been shown to be helpful in eliminating
the so-called high field Q slope, but also reduces the temperature dependent
component of the surface resistance. This treatment is performed right after
the baseline measurement and inside the test cryostat with resistive heaters
strapped to the cavity. In this case, the heaters and temperature sensors are
closer to each other, avoiding a strong temperature swing around the
regulation temperature. The cavity rf space is actively pumped during the bake
and throughout the cryogenic test. During the heat treatment, the surrounding
helium bath space is evacuated to a low mTorr pressure to prevent heat losses
and oxidation of the outer cavity surface. As the bake is done in-situ, the rf
space is not exposed to air after the bake. Witness samples were processed
separately in a sample vacuum oven.

Figure 6:  (a) At 780 MHz and 4 K the mid-T bakes separate further from the
baseline over a larger range of rf fields. (b) At 2 K, the baseline and 120 ∘C
out perform the mid-T bakes.

##  III Results

Figure 7:  (a) At 1166 MHz, the mid-T bakes show a clear improvement in
Q${}_{\textrm{\scriptsize 0}}$ compared to the baseline at 4 K. (b) At 2 K the
400 ∘C bake even out performs the conventional treatments.

### Cavity Performance

Cavity performance measurements were done using standard self-excited loop cw
rf measurements. Both 4.2 K and 2.0 K measurements were performed and the
results of these fixed temperature characterizations are shown in figures 3-7
for all discussed treatments. It is notable that the 120 ∘C/48 hrs bake stands
out as the best performance overall at 4 K and matches the other treatments at
2 K. The mid-T bakes react differently to the change in rf frequency than the
baseline or 120 ∘.

At the lowest frequencies, 220 and 390 MHz (Figs. 3 and 4), the mid-T bakes
match the baseline performance in Q${}_{\textrm{\scriptsize 0}}$ during the 4
K measurements with some marginal changes in quench field. Decreasing the
helium bath temperature to 2 K, the baseline and 120 ∘C treatments show almost
identical performance with G/Q${}_{\textrm{\scriptsize 0}}$ less than 10
n$\Omega$ until high gradient. The two mid-T bakes perform worse both in
Q${}_{\textrm{\scriptsize 0}}$ and quench field, with the 300 ∘C bake
increasing G/Q${}_{\textrm{\scriptsize 0}}$ by about 20 n$\Omega$. The
decrease in Q${}_{\textrm{\scriptsize 0}}$ at these low temperatures and low
frequencies indicates a strong increase in the temperature independent
component R${}_{\textrm{\scriptsize Res}}$ of the surface resistance.

In the next higher frequency modes, 644 and 780 MHz (Figs. 5 and 6), again the
120 ∘C bake is generally the treatment with the highest
Q${}_{\textrm{\scriptsize 0}}$ values, especially at 4 K. Again the mid-T
bakes show a similar performance to the baseline at 4 K, but at low field the
Q${}_{\textrm{\scriptsize 0}}$ is increased. At 780 MHz a larger
Q${}_{\textrm{\scriptsize 0}}$ is measured up to a peak surface field of
around 70 mT. The lower temperature measurement again reveals a poor
performance in the 300 ∘C bake, again indicating a raised
R${}_{\textrm{\scriptsize Res}}$.

At 1166 MHz and 2 K (Fig. 7), the 400 ∘C mid-T bake outperforms the other
treatments. Deconvolution of R${}_{\textrm{\scriptsize s}}$ (Fig. 12) shows
that both mid-T bakes have a much reduced field dependence of
R${}_{\textrm{\scriptsize BCS}}$ but with an increased
R${}_{\textrm{\scriptsize Res}}$. No decrease of either component with
increasing rf amplitude, a signature of mid-T bakes in EP’ed 1.3 GHz
elliptical cavities, is observed.

### Sample Measurements

The witness samples are 12x12 mm${}^{\textrm{\scriptsize 2}}$x1.5 mm flat
coupons made out of high (250+) RRR Nb material. Prior to the witness
treatment they had received a deep etch via BCP and cleaning via an ultrasound
bath. As the 120 ∘C cavity bake was done in the cryostat, the witness sample
processing for this treatment was completed separately from the cavity in a
small sample UHV oven. Therefore, separate samples for the QWR or HWR
treatment only exist for the mid-T bake treatments. These samples were treated
with the cavity together in the induction furnace, but outside the rf volume
due to restricted access into that volume. The same samples are also used in
the $\beta$-SRF beamline at TRIUMF [23] where the magnetic screening
capabilities are being probed.

Figure 8:  EDX data for the baseline, 120 ∘C, and mid-T witness samples. The
400 ∘C bakes shows an enhanced carbon peak with some reduced oxygen for the
HWR. The inset shows the whole energy scan with the Nb peaks.

Energy dispersive x-ray spectroscopy (EDX) is a technique used for an
elemental characterization of the surface of samples. The sample is subjected
to an electron beam. This beam can kick out electrons in the sample atoms. A
higher shell electron drops down, releasing a photon of a specific energy.
This allows an elemental characterization of a sample surface, given that the
photon energies are known. A larger count rate indicates a larger abundance of
this particular element. Fig. 8 shows such a scan for the investigated
samples. Clearly visible are the expected peaks of Nb at around 0.2, 1.9, and
2.2 keV. The differences are in the details of the carbon and oxygen peaks at
0.25 and 0.53 keV, respectively. In the baseline, 120 ∘C, and 300 ∘C samples
any carbon signature is hidden in the tail of the Nb peak, while a noticeable
increase is measurable for the two 400 ∘C samples. Oxygen is slightly
suppressed in the HWR 400 ∘C and QWR 300 ∘C samples. This reduced oxygen peak
is not seen in the QWR 400 ∘C sample. Overall, there seems to be only minor
differences in the carbon and oxygen abundances between the baseline and mid-T
bake samples. Secondary electron microscopy (SEM) of the witness samples show
some contaminates, as is shown in Fig. 9 from the QWR 300 ∘C witness sample.
Similar objects have been seen on the other mid-T baked samples, but not on
the baseline or 120 ∘C sample.

Figure 9:  This SEM image from the QWR 300 ∘C witness sample shows the sampled
area and a contaminating particulate.

The probing depth of EDX is fixed by the energy of the electron beam and depth
information is not obtained. In order to get information about the elemental
composition as a function depth secondary ion mass spectroscopy (SIMS) was
used. In SIMS, the sample is sputtered with a heavy ion beam. For this study
an ION-TOF (GmbH) TOF-SIMS IV instrument supplied a 25 keV Bi${}_{3}^{+}$ beam
pulsed at 10 kHz to bombard the sample surface in an area of about 400x400
$\mu$m2 to generate secondary ions. This releases atoms and molecules from the
sample, which then can be analysed using a time-of-flight mass spectrometer. A
larger sputter time results in a deeper probing of the sample. In Fig. 10 both
the baseline and 400 ∘C QWR witness samples are analysed using this technique.
A larger count rate in the carbon content is noticeable at shallower depths
with the mid-T sample, while the other traces show only minor differences. The
NbO${}_{\textrm{\scriptsize 2}}$, Nb${}_{\textrm{\scriptsize
2}}$O${}_{\textrm{\scriptsize 5}}$, and O${}_{\textrm{\scriptsize 2}}$
signatures in the mid-T bake shift to shallower depths compared to the
baseline sample. Atomic oxygen shows similar count rates between the two
samples. At larger sample depths, all traces converge to similar values,
showing that the heat treatment changes only the surface composition.

A similar shift in the Nb${}_{\textrm{\scriptsize
2}}$O${}_{\textrm{\scriptsize 5}}$ layer has been observed on samples after a
mid-T bake [24]. Increased carbon impurities had also been observed [11] in
mid-T samples and associated with weaker rf performance.

Figure 10: SIMS data of the baseline (dashed lines) and 400 ∘C (solid lines)
QWR witness sample for a number of elements of interest. A notable change is
in the carbon trace with a larger count rate in the mid-T sample.

##  IV Discussion

The rf performance of both cavities shows clear changes between the various
heat treatments, frequencies, and temperatures. To analyse the performance,
the surface resistance is deconvoluted into its components
R${}_{\textrm{\scriptsize BCS}}$ and R${}_{\textrm{\scriptsize Res}}$ by
fitting R${}_{\textrm{\scriptsize s}}$(T) data collected during the cooldown
from 4.2 K to 2 K at fixed rf peak fields B${}_{\textrm{\scriptsize p}}$ from
10 up to 110 mT to Eq. 2. Results of these fits are shown in figures 11 and 12
for the 220 and 1166 MHz modes. The baseline and 120 ∘C treatments show in all
frequencies and temperatures a strong rf field dependence of
R${}_{\textrm{\scriptsize BCS}}$, with the 120∘C bake producing significantly
lower resistances. The residual resistance R${}_{\textrm{\scriptsize Res}}$ is
marginally increased as is commonly observed for this treatment. The two
investigated mid-T bakes mostly follow similar trends except at low
temperatures and high frequency (Fig. 12b). There, R${}_{\textrm{\scriptsize
BCS}}$ has almost no field dependence and is lower than for the 120∘C bake. A
higher Q${}_{\textrm{\scriptsize 0}}$ compared to the 120 ∘C bake is prevented
by a larger R${}_{\textrm{\scriptsize Res}}$, which could be a consequence of
the higher carbon content as measured by SIMS.

To quantify the dependence of R${}_{\textrm{\scriptsize BCS}}$ on rf field
amplitude B${}_{\textrm{\scriptsize p}}$, the data is fitted with a quadratic
function according to

$\displaystyle
R_{BCS}=R_{0}\left(1+\gamma\left(\frac{B_{p}}{B_{0}}\right)^{2}\right)$ (6)

with $R_{0}$ as zero-field resistance, $B_{0}$ as a normalizing constant
(chosen to be 100 mT) and $\gamma$ as field dependency. The quadratic form
originates from thermal feedback models [25] and other models [26] and is here
adopted to quantify the changes in field dependence and zero-field resistance
of R${}_{\textrm{\scriptsize BCS}}$. Fig. 13a shows R${}_{\textrm{\scriptsize
0}}$ at 2 K resolved as function of rf frequency $f$ for the different
investigated treatments. For all measured treatments and cavity modes, the
baseline treatment results in the highest R${}_{\textrm{\scriptsize 0}}$ while
the 120 ∘C bake produces the lowest zero-field resistance. The mid-T bakes
fall in between.

The baseline and 120 ∘C treatments show a steady increasing
R${}_{\textrm{\scriptsize BCS}}$ with frequency as would be expected from eq.
2, with the low temperature bake showing a lower slope, thus a reduced
frequency dependence. It is noteworthy that the 400 ∘C treatment shows a lower
resistance at 1166 MHz than at 780 MHz, which is an indication that this heat
treatment is only effective at high rf frequencies above 1 GHz. The 300 ∘C
bake shows a frequency dependence of R${}_{\textrm{\scriptsize 0}}$ comparable
to the 120 ∘C bake.

The rf amplitude dependence $\gamma$ for 2 K shows no clear correlation with
treatment or rf frequency as can be seen in Fig. 13(b). For the baseline and
mid-T treatments, $\gamma$ stays roughly constant and at similar values for
all but the lowest frequency mode. The 220 MHz mode shows the highest field
dependence for these three treatments. At this frequency, $\gamma$ is
significantly higher for the baseline and 400 ∘C bake compared to the other
resonant modes. The 120 ∘C bake on the other hand has its highest $\gamma$ at
1166 MHz, while $\gamma$ for the other frequencies stays constant. This
together with the changes in R${}_{\textrm{\scriptsize 0}}$ by the treatments,
shows that the differences in R${}_{\textrm{\scriptsize BCS}}$ between
treatments are mainly caused by the zero-field resistance
R${}_{\textrm{\scriptsize 0}}$.

Figure 11: The 120 ∘C bake outperforms the other treatments at low rf
frequency, 220 MHz in this case, due to a reduced R${}_{\textrm{\scriptsize
BCS}}$ ((a) 4K and (b) 2K), while increasing R${}_{\textrm{\scriptsize Res}}$
(c).

Using all five available modes from the two cavities, the frequency dependence
of R${}_{\textrm{\scriptsize BCS}}$ is extracted to test against the
theoretical prediction given by Eq. 2 and codes like SRIMP which is based on
[17]. A power law of the form

$\displaystyle R_{BCS}(\omega)=A\omega^{x}$ (7)

can be used to model the relationship with $A$ and $x$ as fit parameters.
Based on SHRIMP, the theoretical prediction for $x$ is around 1.87. $x$ as a
function of field amplitude is shown in figure 14 for 4 K and 2 K. The
baseline treatment matches the prediction closest for a wide range of rf
fields and at both temperatures. The other treatments deviate from the
prediction as rf amplitude increases. At 4 K, the baseline, 120 ∘C and 400 ∘C
all follow a similar downward trend, with the 300 ∘C treatment resulting in
significantly lower values than the other treatments. At 2 K, both mid-T
treatments show a reduced frequency dependent R${}_{\textrm{\scriptsize BCS}}$
compared to the baseline, 120 ∘C bake or the theoretical prediction. This is
an indicator of a changed layer composition in the oxygen diffused cavities,
so that the resistance deviates from the typical assumed Mattis-Bardeen
theory. It should be noted that the theory, and the resulting model for the
surface resistance, is a low field approximation, and may not be valid at
higher rf fields.

The sample analysis done with SIMS shows a modification of the oxide layer due
to the mid-T bake, which is consistent with previous reports [9, 11, 24]. The
significant CO${}_{\textrm{\scriptsize 2}}$ partial pressure during the ramp
up of the heat treatment (Fig. 2) and the strong carbon concentration in the
SIMS measurement indicate that the cavity surfaces could have been
contaminated at that time. It is speculated that this contamination causes the
high R${}_{\textrm{\scriptsize Res}}$ in the mid-T baked cavity measurements.

Figure 12:  The mid-T bakes outperform in R${}_{\textrm{\scriptsize BCS}}$ the
baseline and 120 ∘C bake at high rf frequency, 1166 MHz shown here, and low
temperature (b) but a large R${}_{\textrm{\scriptsize Res}}$ (c) prevents a
gain in Q${}_{\textrm{\scriptsize 0}}$.

The changed frequency response of R${}_{\textrm{\scriptsize BCS}}$ with the
mid-T bakes shows that the treatments have an effect on the cavity performance
and that these treatments are strongly rf frequency dependent. This effect and
the lack of a decreasing R${}_{\textrm{\scriptsize BCS}}$ is similar to the
frequency response of N-doped cavities as EP’ed 650 MHz cavities do not show a
decreasing R${}_{\textrm{\scriptsize BCS}}$ [13]. A similar observation was
made in a 650 MHz EP’ed cavity with a 300 ∘C furnace bake [14]. Only high rf
frequencies above 1 GHz seem to be able to produce this characteristic. In
addition, to get the highest Q${}_{\textrm{\scriptsize 0}}$ out of N-doping,
EP is required over BCP [12]. The grain structure seems to play a role as fine
grain cavities, which are most commonly used, treated with BCP did not show
the otherwise typical rise in Q${}_{\textrm{\scriptsize 0}}$ even at 1.3 GHz.
Mid-T furnace baked cavities at 1.3 GHz show a similar requirement [11].

Figure 13:  (a) R${}_{\textrm{\scriptsize 0}}$ at 2 K increases as the
frequency increases. The 120 ∘C bake has a significantly reduced
R${}_{\textrm{\scriptsize 0}}$ compared to the other treatments. Of note is
also that the 400 ∘C bake shows a reduced R${}_{\textrm{\scriptsize 0}}$ at
high frequency. (b) The field dependence as determined by $\gamma$ on the
other hand has no clear correlation with frequency.

This all leads to the conclusion that high frequency above 1 GHz and electro-
polishing are required to measure decreasing surface resistances in oxygen
diffused cavities.

##  V Summary

In this paper, the performance of two coaxial multi-mode cavities after four
different treatments is presented giving a unique insight into the frequency
dependence of R${}_{\textrm{\scriptsize BCS}}$. The investigated baseline and
120 ∘C bakes are established treatments, while mid-T bakes are novel and have
not been reported on for low frequency coaxial type cavities before. Mid-T
bake recipes have been shown to produce cavity performance similar to
N${}_{\textrm{\scriptsize 2}}$ doped cavities. Two different mid-T treatment
temperatures, 400 and 300 ∘C, were applied to the QWR and HWR multi-mode
cavities and their performance measured at 4.2 and 2.0 K. In addition,
Q${}_{\textrm{\scriptsize 0}}$ measurements during the cooldown were done to
extract the temperature dependence of the surface resistance. Measurements
were done in five resonant modes, spanning from 200 MHz to 1.2 GHz, covering a
wide range of relevant frequencies for SRF particle accelerators.

In the low frequency modes at 220 and 390 MHz, the mid-T bake treatments at
400 and 300 ∘C did not improve the performance over either the baseline or the
120 ∘C bake. At 4 K, the 120 ∘C bake produces the highest
Q${}_{\textrm{\scriptsize 0}}$ in all cases and therefore the lowest
R${}_{\textrm{\scriptsize BCS}}$. As rf frequency increases, the mid-T bakes
separate from the baseline performance and perform better. At 1.2 GHz and 2 K,
R${}_{\textrm{\scriptsize BCS}}$ is lower than what was measured for the 120
∘C bake, but no decreasing resistance with rf field could be observed as it
has been in elliptical EP’ed cavities. There is a clear difference in the
frequency dependence between the established treatments and the mid-T bakes,
especially at higher rf fields as figure 14 shows. All this leads to the
conclusion that there is a strong impact of the rf frequency on the
performance of cavities with a mid-T bake.

Figure 14:  The frequency dependence exponent $x$ (eq. 7) at (a) 4 K and (b) 2
K as a function of rf field amplitude. At 4 K and higher rf fields the 300 ∘C
bake shows a much lower frequency dependence than the other treatments, while
at the lower temperature the trends separate between the established
treatments and the newer mid-T bakes.

Witness sample measurements using EDX and SIMS show some changes in the
elemental composition and oxide layers. EDX shows that some of the mid-T
samples have a heightened carbon content on the surface. This is confirmed
with the SIMS measurement of a 400 ∘C mid-T witness sample. A source of the
carbon could be the UHV furnace used for the mid-T bakes as the RGA spectrum
during the heat treatment shows a significant CO${}_{\textrm{\scriptsize 2}}$
peak during the temperature ramp up as shown in figure 2. In addition, changes
to the oxide layer were measured which are consistent with the oxide diffusion
treatment. The lack of rising Q${}_{\textrm{\scriptsize 0}}$ in these BCP’ed
cavities leads to the hypothesis that electro-polishing, and the resulting
smoother surface, is needed to reach the full potential of this treatment.
This might be a similar phenomenon to N${}_{\textrm{\scriptsize 2}}$ doping
where BCP’ed cavities do not show the increase in Q${}_{\textrm{\scriptsize
0}}$.

## Acknowledgements

The authors would like to thank the technical members of the SRF group, namely
Devon Lang, James Keir, Bhalwinder Waraich, Ruminder Sekhon, and cryogenics
group, namely David Kishi, Johnson Cheung, Rowan Bjarnson, at TRIUMF who
helped tremendously with performing the cavity preparations and tests, as well
as the UWO team for the SIMS measurements. Funding is provided by Natural
Sciences and Engineering Research Council of Canada (NSERC) grants
SAPPJ-2020-00030 (holder T. Junginger) and SAPIN-2019-00040, SAPPJ-2017-00035,
and SAPEQ-2015-00005 (holder R.E. Laxdal).

## References

  * Dilling et al. [2014] Dilling J, Krücken R, Ball G.  Isac overview.  Hyperfine Interactions 225 (2014).  doi:10.1007/s10751-013-0877-7.
  * Grassellino et al. [2017] Grassellino A, Romanenko A, Trenikhina Y, Checchin M, Martinello M, Melnychuk OS, et al.  Unprecedented quality factors at accelerating gradients up to 45 MVm-1in niobium superconducting resonators via low temperature nitrogen infusion.  Superconductor Science and Technology 30 (2017) 094004.  doi:10.1088/1361-6668/aa7afe.
  * Grassellino et al. [2013] Grassellino A, Romanenko A, Sergatskov D, Melnychuk O, Trenikhina Y, Crawford A, et al.  Nitrogen and argon doping of niobium for superconducting radio frequency cavities: a pathway to highly efficient accelerating structures.  Superconductor Science and Technology 26 (2013) 102001.  doi:10.1088/0953-2048/26/10/102001.
  * Gonnella et al. [2019] Gonnella D, Aderhold S, Bafia D, Burrill A, Checchin M, Ge M, et al.  The LCLS-II HE High Q and Gradient R&D Program.  Proc. SRF’19 (JACoW Publishing, Geneva, Switzerland) (2019).  doi:10.18429/JACoW-SRF2019-MOP045.  Https://doi.org/10.18429/JACoW-SRF2019-MOP045.
  * He et al. [2021] He F, Pan W, Sha P, Zhai J, Mi Z, Dai X, et al.  Medium-temperature furnace baking of 1.3 ghz 9-cell superconducting cavities at ihep.  Superconductor Science and Technology 34 (2021) 095005.  doi:10.1088/1361-6668/ac1657.
  * Ito et al. [2021] Ito H, Araki H, Takahashi K, Umemori K.  Influence of furnace baking on Q–E behavior of superconducting accelerating cavities.  Progress of Theoretical and Experimental Physics 2021 (2021).  doi:10.1093/ptep/ptab056.  071G01.
  * Steder et al. [2022] Steder L, Bate C, Remde H, Reschke D, Schaffran J, Trelle L, et al.  Medium Temperature Treatments of Superconducting Radio Frequency Cavities at DESY.  Proc. 31st International Linear Accelerator Conference (LINAC’22) (JACoW Publishing, Geneva, Switzerland) (2022), no. 31 in International Linear Accelerator Conference, 840–843.  doi:10.18429/JACoW-LINAC2022-THPOGE22.
  * Dhakal [2020] Dhakal P.  Nitrogen doping and infusion in srf cavities: A review.  Physics Open 5 (2020) 100034.  doi:https://doi.org/10.1016/j.physo.2020.100034.
  * Lechner et al. [2021] Lechner EM, Angle JW, Stevie FA, Kelley MJ, Reece CE, Palczewski AD.  RF surface resistance tuning of superconducting niobium via thermal diffusion of native oxide.  Applied Physics Letters 119 (2021).  doi:10.1063/5.0059464.  082601\.
  * Semione et al. [2019] Semione GDL, Pandey AD, Tober S, Pfrommer J, Poulain A, Drnec J, et al.  Niobium near-surface composition during nitrogen infusion relevant for superconducting radio-frequency cavities.  Phys. Rev. Accel. Beams 22 (2019) 103102.  doi:10.1103/PhysRevAccelBeams.22.103102.
  * Yang et al. [2022] Yang Z, Hao J, Quan S, Lin L, Wang F, Jiao F, et al.  Surface resistance effects of medium temperature baking of buffered chemical polished 1.3 ghz nine-cell large-grain cavities.  Superconductor Science and Technology 36 (2022) 015001.  doi:10.1088/1361-6668/aca12a.
  * Dhakal et al. [2019] Dhakal P, Ciovati G, Pudasaini U, Chetri S, Balachandran S, Lee PJ.  Surface characterization of nitrogen-doped high purity niobium coupons compared with superconducting rf cavity performance.  Phys. Rev. Accel. Beams 22 (2019) 122002.  doi:10.1103/PhysRevAccelBeams.22.122002.
  * Martinello et al. [2021] Martinello M, Bice DJ, Boffo C, Chandrasekeran SK, Eremeev GV, Furuta F, et al.  Q-factor optimization for high-beta 650 MHz cavities for PIP-II.  Journal of Applied Physics 130 (2021).  doi:10.1063/5.0068531.  174501\.
  * Sha et al. [2022] Sha P, Pan W, Zhai J, Mi Z, Jin S, Liu B, et al.  Quality factor enhancement of 650 mhz superconducting radio-frequency cavity for cepc.  Applied Sciences 12 (2022).  doi:10.3390/app12020546.
  * Kolb et al. [2020] Kolb P, Yao Z, Junginger T, Dury B, Fothergill A, Vanderbanck M, et al.  Coaxial multimode cavities for fundamental superconducting rf research in an unprecedented parameter space.  Phys. Rev. Accel. Beams 23 (2020) 122001.  doi:10.1103/PhysRevAccelBeams.23.122001.
  * Gurevich [2017] Gurevich A.  Theory of RF superconductivity for resonant cavities.  Superconductor Science and Technology 30 (2017) 034004.  doi:10.1088/1361-6668/30/3/034004.
  * Halbritter [1970] Halbritter J.  Fortan-program for the computation of teh surface impedance of superconductors.  KFK-external report 3/70-6 (1970).
  * Delayen et al. [2018] Delayen JR, Park H, De Silva SU, Ciovati G, Li Z.  Determination of the magnetic field dependence of the surface resistance of superconductors from cavity tests.  Phys. Rev. Accel. Beams 21 (2018) 122001.  doi:10.1103/PhysRevAccelBeams.21.122001.
  * Palmieri et al. [2014] Palmieri V, Rossi AA, Stark SY, Vaglio R.  Evidence for thermal boundary resistance effects on superconducting radiofrequency cavity performances.  Superconductor Science and Technology 27 (2014) 085004.  doi:10.1088/0953-2048/27/8/085004.
  * Dhakal et al. [2012] Dhakal P, Ciovati G, Rigby W, Wallace J, Myneni GR.  Design and performance of a new induction furnace for heat treatment of superconducting radiofrequency niobium cavities.  Review of Scientific Instruments 83 (2012) 065105.  doi:10.1063/1.4725589.
  * Angle et al. [2023] Angle JW, Lechner EM, Reece CE, Stevie FA, Kelley MJ.  Analysis of furnace contamination on superconducting radio frequency niobium using secondary-ion mass spectrometry.  Journal of Vacuum Science & Technology B 41 (2023).  doi:10.1116/6.0002624.  034005\.
  * Ciovati [2004] Ciovati G.  Effect of low-temperature baking on the radio-frequency properties of niobium superconducting cavities for particle accelerators.  Journal of Applied Physics 96 (2004) 1591–1600.  doi:10.1063/1.1767295.
  * Thoeng et al. [2023] Thoeng E, McFadden RML, Saminathan S, Morris GD, Kolb P, Matheson B, et al.  A new high parallel-field spectrometer at TRIUMF’s $\beta$-NMR facility.  Review of Scientific Instruments 94 (2023).  doi:10.1063/5.0137368.  023305\.
  * Posen et al. [2020] Posen S, Romanenko A, Grassellino A, Melnychuk O, Sergatskov D.  Ultralow surface resistance via vacuum heat treatment of superconducting radio-frequency cavities.  Phys. Rev. Appl. 13 (2020) 014024.  doi:10.1103/PhysRevApplied.13.014024.
  * Gurevich [2006] Gurevich A.  Thermal rf breakdown of superconducting cavities.  ICFA Beam Dyn. Newslett 39 (2006) 34–44.
  * Weingarten [2011] Weingarten W.  Field-dependent surface resistance for superconducting niobium accelerating cavities.  Phys. Rev. ST Accel. Beams 14 (2011) 101002.  doi:10.1103/PhysRevSTAB.14.101002.