11institutetext: Departament de Física, Universitat de les Illes Balears,
E-07122, Palma de Mallorca, Spain
22institutetext: Insitut d’Aplicacions Computacionals de Codi Comunitari
(IAC3), Universitat de les Illes Balears, E-07122, Palma de Mallorca, Spain

# Self-consistent equilibrium models of prominence thin threads heated by
Alfvén waves propagating from the photosphere

Llorenç Melis 1122 Roberto Soler 1122 Jaume Terradas 1122

The fine structure of solar prominences is made by thin threads that outline
the magnetic field lines. Observations show that transverse waves of Alfvénic
nature are ubiquitous in prominence threads. These waves are driven at the
photosphere and propagate to prominences suspended in the corona. Heating due
to Alfvén wave dissipation could be a relevant mechanism in the cool and
partially ionised prominence plasma. Here, we explore the construction of 1D
equilibrium models of prominence thin threads that satisfy energy balance
between radiative losses, thermal conduction, and Alfvén wave heating. We
assume the presence of a broadband driver at the photosphere that launches
Alfvén waves towards the prominence. An iterative method is implemented, in
which the energy balance equation and the Alfvén wave equation are
consecutively solved. From the energy balance equation and considering no wave
heating initially, we compute the equilibrium profiles along the thread of the
temperature, density, ionisation fraction, etc. On these equilibrium profiles,
we use the Alfvén wave equation to compute the wave heating rate, which is
then put back in the energy balance equation to obtain new equilibrium
profiles, and so on. The process is repeated until convergence to a self-
consistent thread model heated by Alfvén waves is achieved. We have obtained
equilibrium models composed of a cold and dense thread, a extremely thin PCTR,
and an extended coronal region. The length of the cold thread decreases with
the temperature at the prominence core and increases with the Alfvén wave
energy flux injected at the photosphere. However, equilibrium models are not
possible for sufficiently large wave energy fluxes when the wave heating rate
inside the cold thread becomes larger than radiative losses. The maximum value
of the wave energy flux that allows an equilibrium depends on the prominence
core temperature. This constrains the existence of thread equilibria in
realistic conditions.

###### Key Words.:

magnetohydrodynamics (MHD) – Sun:atmosphere – Sun:corona –
Sun:filaments,prominences – Sun:oscillations – waves

##  1 Introduction

Solar prominences consist of masses of relatively cool and dense plasma
suspended in the solar corona, whose physical properties are similar to those
in the chromosphere (see, e.g., Vial & Engvold, 2015). High-resolution
observations have shown that prominences are composed by a myriad of thin
threads, which seem to outline particular field lines of the general magnetic
structure of prominences (see, e.g., Lin, 2011; Martin, 2015). The mechanical
equilibrium of prominences suspended above the photosphere can be established
when the upward force provided by the magnetic field compensates gravity (see,
e.g., Parenti, 2014; Gilbert, 2015; Heinzel, 2015). On the other hand, the
energy balance in solar prominences is a problem not so well understood
(Gilbert, 2015). A detailed knowledge of the heating and cooling processes
operating in prominences is needed to explain the temperatures at prominence
cores, which are estimated to be in the range 7,000–9,000 K (see Parenti,
2014). According to radiative-equilibrium models, heating due to radiative
illumination is believed to be the main heating source (e.g., Heasley &
Mihalas, 1976; Anzer & Heinzel, 1999; Heinzel et al., 2010; Heinzel & Anzer,
2012). Using a slab model, Heinzel & Anzer (2012) computed radiative-
equilibrium temperatures within the expected range of values for a plasma
composed of hydrogen alone. However, if CaII losses are added, much lower
radiative-equilibrium temperatures are obtained and, for certain conditions,
they can be as low as 4,400 K (see details in Heinzel & Anzer, 2012). To the
best of our knowledge, the effect of additional coolants has not been
explored, but it may be that even lower radiative-equilibrium temperatures
would be obtained if additional coolant species are added. Therefore, although
radiative heating would be dominant in prominences, we cannot presently
discard that additional sources of heating are also playing a role. The
purpose of the present paper is to explore the role of Alfvén wave dissipation
as a potential heating mechanism in solar prominences that could complement
radiative heating.

Mechanical and thermal equilibrium models of prominence thin threads have been
studied before. For instance, Degenhardt & Deinzer (1993) studied the
equilibrium between cooling and an imposed heating without considering the
effect of thermal conduction. More detailed equilibrium models of threads were
constructed by Terradas et al. (2021), who studied the balance between
radiative cooling, thermal conduction, and heating. However, the heating
considered by Terradas et al. (2021) was arbitrarily imposed rather than being
the result of an actual physical process. The present work follows the method
of Terradas et al. (2021) to construct equilibrium models but using a
consistently computed heating from the dissipation of Alfvén waves.

Observations have shown that transverse waves are ubiquitous in fine
structures of solar prominences (see the review by Arregui et al., 2018).
These waves are interpreted as magnetohydrodynamic waves of Alfvénic nature
(see, e.g., Lin et al., 2009; Ballester, 2015). There are strong indications
that the transverse waves observed in prominences are driven at the
photosphere, where the prominence magnetic field is anchored (Hillier et al.,
2013). It has been shown that Alfvénic waves can travel from the photosphere
to the coronal structures, like prominences or coronal loops, transporting a
significant amount of energy (see Soler et al., 2019, 2021). The existence of
an efficient dissipation mechanism in the prominence plasma could provide a
way to thermalise the wave energy. The fact that the prominence plasma is only
partially ionised introduces important dissipation mechanisms for Alfvén
waves, namely ambipolar diffusion due to ion-neutral collisions and enhanced
Ohm’s diffusion due to electron-neutral collisions (see Ballester et al.,
2018). Thus, heating due to Alfvénic waves arises as a possible heating
mechanism in solar prominences whose importance needs to be explored.

There are some previous works that have explored the role of Alfvén wave
heating in prominences. Using monochromatic waves and considering ion-neutral
collisions, Pécseli & Engvold (2000) concluded that Alfvén wave heating may
only compensate for a tiny part of radiative losses. However, these authors
did not consider a more realistic broadband spectrum of waves, whose presence
is observationally confirmed (Hillier et al., 2013). Conversely, Parenti &
Vial (2007) estimated that the heating produced by the dissipation of
unresolved non-thermal motions in prominences, which were attributed to Alfvén
waves, could indeed compensate for a large fraction of radiative losses.
Later, Soler et al. (2016) used an idealised slab model with a straight
magnetic field transverse to the slab axis. The slab was filled with a
homogeneous and partially ionised prominence plasma and was embedded in a
fully ionised corona. No fine structure was considered. Alfvén waves were
launched towards the prominence slab and heating due to ion-neutral collisions
was computed. Soler et al. (2016) estimated that the heating could account for
about 10 $\%$ of radiative losses if a spectrum of wave periods between 0.1 to
100 s is considered. Although the model was simple, the results of Soler et
al. (2016) evidenced the potential of wave heating in solar prominences.

A further study of wave heating in solar prominences was done in Melis et al.
(2021), where the model implemented was more elaborated. An important
difference with the previous work by Soler et al. (2016), is that Melis et al.
(2021) considered a model for a thin thread including longitudinal non-
uniformity and the anchoring of the magnetic field at the ends of the model
representing the base of the corona. In the central part of the thread, where
the plasma is densest and coolest, partial ionization and the roles of Ohm’s
and ambipolar diffusion were included. In the model of Melis et al. (2021),
Alfvén waves were driven at one end on the thread to mimic photospheric
excitation. Their results showed that wave heating could compensate for
radiative losses in the cool part of the thread, where the plasma was
partially ionised, although the energy balance was not consistently satisfied
because the density and temperature profiles were assumed ad hoc and the back
reaction of wave heating upon those profiles was not considered.

The aim of the present work is to go further and continue the study of Alfvén
wave heating in prominence thin threads. Here the goal is to investigate the
effect of wave heating in the equilibrium of threads. To this end, we shall
combine the method to compute the wave heating rate of Melis et al. (2021)
with the method to compute equilibrium models of Terradas et al. (2021). In
the presence of wave heating, threads should tend to an equilibrium
configuration where wave heating consistently affects the profiles of density,
temperature, and other variables along the thread. To study such influence and
to investigate under what conditions an equilibrium is possible, we will
implement a self-consistent approach loosely inspired by that used in Ofman et
al. (1998) for heating by resonant absorption in coronal loops. Starting from
a thread model that includes no wave heating, the wave heating rate will be
computed following Melis et al. (2021). Then, following Terradas et al.
(2021), the energy balance equation will be solved including wave heating, and
a new thread model will be obtained. The wave heating rate will be re-
calculated in the new thread model and energy balance will be imposed again.
This, in turn, will lead to a different thread model, and so on. This scheme
will be run iteratively until convergence to a self-consistent model is
reached.

This paper is structured as follows. Section 2 contains an explanation of the
thread model, the basic equations, and the implementation of the self-
consistent method. In Section 3, we present and discuss the results. Finally,
some conclusions are given in Section 4.

##  2 Method

###  2.1 Thread models satisfying energy balance

As in Melis et al. (2021), we use a 1D prominence thread model composed of a
single magnetic field line of length $L=10^{8}$ m. Curvature and gravity are
not considered, hence the magnetic field line is straight with a uniform value
of $B=10$ G. We align the magnetic field with the $z$-axis for convenience.
The centre of the thread is located at $z=0$, whereas its ends are located at
$z=\pm L/2$. For simplicity, we assume a uniform gas pressure of $p=5\times
10^{-3}$ Pa. A sketch of the thread model is shown in Figure 1.

Figure 1: Sketch of the 1D thread model. The grey areas denote the
photosphere, where the thread ends are anchored. The red color gradient aims
to represent the density distribution along the thread.

The density, $\rho$, and temperature, $T$, profiles are expected to be non-
uniform along the thread, but $p$ is assumed to be uniform in order to
preserve the hydrostatic equilibrium. Density and temperature are related with
the pressure through the ideal gas law, which for a hydrogen plasma is

$p=(1+\xi_{i})\rho RT,$ (1)

where $\xi_{i}$ is the ionisation fraction defined as the ratio between ion
density with total density and $R$ is the ideal gas constant. In the cool
parts of the thread, we expect the plasma to be partially ionised and so
$\xi_{i}<1$. In order to realistically obtain the ionisation fraction in a
solar prominence thread, we would need to solve the full radiative transfer
equations including incident radiation. This is far beyond the aims of this
work. Instead, we approximate the ionization degree using the tabulated values
given in Heinzel et al. (2015), which are based on non-LTE computations for 1D
prominence slabs. In Table 1 we indicate the values of the ionisation fraction
as a function of the temperature for the pressure value considered here. The
data of Table 1 is interpolated for intermediate temperatures and $\xi_{i}=1$
is imposed for $T\geq$ 20,000 K. If either $T$ or $\rho$ is known, Equation
(1) can be used to compute the other variable. Hence, we shall focus on the
computation of $T$.

$T$ | 6000 | 8000 | 10000 | 12000 | 14000
---|---|---|---|---|---
$\xi_{i}$ | 0.41 | 0.52 | 0.68 | 0.81 | 0.89
Table 1: Values of the ionization fraction, $\xi_{i}$, for various values of
the temperature (in K) for a gas pressure of $5\times 10^{-3}$ Pa and an
altitude above the photosphere of 20,000 km. Adapted from Heinzel et al.
(2015)

In order to obtain the temperature profile, we solve the energy balance
equation given by

$\langle Q\rangle-\nabla\cdot\vec{q}-\mathcal{L}+C=0,$ (2)

where $\langle Q\rangle$ denotes the heating rate (this will be the wave
heating, as explained later), $\vec{q}=-\kappa_{\parallel}\nabla T$ is the
heat flux due to thermal conduction with $\kappa_{\parallel}$ the parallel
thermal conductivity, $\mathcal{L}$ represents the energy losses caused by
radiative cooling, and $C$ may represent other heating sources like, e.g.,
heating due to viscous dissipation or radiative heating. Alfvén wave
dissipation due to viscosity is less efficient than ambipolar dissipation in
prominences, while a consistent treatment of radiative heating is well beyond
our aims. Therefore, these additional sources are ignored here.

In a partially ionised plasma, the parallel component of the thermal
conductivity with respect to the magnetic field can be approximated by
$\kappa_{\parallel}\approx\kappa_{e}+\kappa_{n}$, where $\kappa_{e}$ and
$\kappa_{n}$ are the contributions of electrons and neutrals, respectively,
which are given by

$\displaystyle\kappa_{e}$ $\displaystyle=$ $\displaystyle
3.2\frac{n_{e}^{2}k_{B}^{2}T}{\alpha_{e}+\alpha_{ee}},$ (3)
$\displaystyle\kappa_{n}$ $\displaystyle=$
$\displaystyle\frac{5}{3}\frac{n_{n}^{2}k_{B}^{2}T}{\alpha_{n}+\alpha_{nn}},$
(4)

where $k_{B}$ is the Boltzmann constant, $n_{e}$ and $n_{n}$ are the number
density of electrons and neutrals, respectively, $\alpha_{e}$ and $\alpha_{n}$
are the corresponding electron and neutral total friction coefficients, and
$\alpha_{ee}$ and $\alpha_{nn}$ are the electron and neutral friction
coefficients accounting for self-collisions. The total friction coefficients
are computed as

$\alpha_{\beta}=\sum_{\beta\neq\beta^{\prime}}\alpha_{\beta\beta^{\prime}},$
(5)

where $\beta$ and $\beta^{\prime}$ denote the different species that collide,
that could be electrons (e), neutrals (n) or ions (i), and
$\alpha_{\beta\beta^{\prime}}=\alpha_{\beta^{\prime}\beta}$ are the symmetric
friction coefficients between individual particles. If the collision is
between charged particles, the friction coefficient is

$\alpha_{\beta\beta^{\prime}}=\frac{n_{\beta}n_{\beta^{\prime}}e^{4}\ln{\Lambda_{\beta\beta^{\prime}}}}{6\pi\sqrt{2\pi}\epsilon_{0}^{2}m_{\beta\beta^{\prime}}\left(k_{B}T/m_{\beta\beta^{\prime}}\right)^{3/2}},$
(6)

where
$m_{\beta\beta^{\prime}}=m_{\beta}m_{\beta^{\prime}}/\left(m_{\beta}+m_{\beta^{\prime}}\right)$
is the reduced mass, $e$ is the electron charge, $\epsilon_{0}$ is the
electrical permittivity, $n_{\beta}$ is the number density, and
$\ln\Lambda_{\beta\beta^{\prime}}$ is the Coulomb’s logarithm, expressed as

$\ln{\Lambda_{\beta\beta^{\prime}}}=\ln{\left(\frac{24\pi\epsilon_{0}^{3/2}k_{B}^{3/2}T^{3/2}}{e^{3}\sqrt{n_{\beta}+n_{\beta^{\prime}}}}\right)}.$
(7)

If at least one particle is neutral, the friction coefficient is

$\alpha_{\beta n}=n_{\beta}n_{n}m_{\beta n}\left[\frac{8k_{B}T}{\pi m_{\beta
n}}\right]^{1/2}\sigma_{\beta n},$ (8)

where $\sigma_{\beta n}$ is the collision cross-section. See details in Melis
et al. (2021).

The radiative cooling rate is computed using the function provided in Athay
(1986), which is frequently used in the literature to approximate radiative
losses in cool prominence plasmas (see, e.g., Mok et al., 1990). We use
Athay’s function because it provides lower radiative losses than other cooling
tables based on the optically thin approximation that may overestimate the
actual losses for cool prominence temperatures. However, we note that other
approaches are possible (see Dalgarno & McCray, 1972), which more recent
cooling functions implement (e.g., Schure et al., 2009; Hermans & Keppens,
2021; Brughmans et al., 2022). The adopted radiative function is written down
as

$\mathcal{L}(\rho,T)=f_{p}(T)\frac{\rho^{2}T^{2}}{m_{p}^{2}},$ (9)

where $m_{p}$ is the proton mass and $f_{p}$ is an analytical function of the
temperature that can be cast in MKS units as

$\displaystyle f_{p}(T)$ $\displaystyle=$ $\displaystyle
10^{-35}T^{-2}\Bigl{\\{}0.4\exp\left[-30\left(\log_{10}T-4.6\right)^{2}\right]$
(10) $\displaystyle+4\exp\left[-20\left(\log_{10}T-4.9\right)^{2}\right]$
$\displaystyle+4.5\exp\left[-16\left(\log_{10}T-5.35\right)^{2}\right]$
$\displaystyle+2\exp\left[-4\left(\log_{10}T-6.1\right)^{2}\right]\Bigr{\\}}.$

We note that because of the variation along the thread of the temperature,
density, and ionization degree, both the thermal conductivity,
$\kappa_{\parallel}$, and the radiative cooling rate, $\mathcal{L}$, are
functions of $z$. For a given profile of the heating rate, $\langle Q\rangle$,
Equation (2) results in a complicated second order ordinary differential
equation for the temperature profile, $T$, that needs to be solved
numerically. To do so, we consider the same boundary conditions as in Terradas
et al. (2021). We prescribe the temperature at the thread centre, $T_{0}$, and
impose that this must correspond to the temperature minimum, so that the
boundary conditions are

$\displaystyle T=T_{0},\qquad\textrm{at}\qquad z=0,$ (11)
$\displaystyle\frac{\partial T}{\partial z}=0,\qquad\textrm{at}\qquad z=0.$
(12)

Terradas et al. (2021) showed that the following relation must be satisfied at
the thread centre,

$\frac{\partial^{2}T}{\partial z^{2}}=\frac{\mathcal{L}-\langle
Q\rangle}{\kappa_{\parallel}}>0,\qquad\textrm{at}\qquad z=0.$ (13)

For a cold thread, the temperature at the centre must be a minimum. According
to Equation (13), this is satisfied only when $\mathcal{L}>\langle Q\rangle$
at $z=0$. This gives us a restriction on the value of the heating rate that
will be explored further.

The numerical integration of Equation (2) is performed in two separate stages,
first from $z=0$ to $z=L/2$ and later from $z=0$ to $z=-L/2$. Then, the two
solutions are joined together to construct the whole profile. This is done to
correctly account for asymmetries in the heating rate, $\langle Q\rangle$. The
integration is performed in _Wolfram Mathematica_ with the routine `NDSolve`,
which automatically adapts the step size to minimise numerical errors.

###  2.2 Alfvén wave heating

As in Melis et al. (2021), the heating rate $\langle Q\rangle$ is computed
from the dissipation of Alfvén waves. In the single-fluid magnetohydrodynamic
approximation and assuming linear waves, an equation that governs the
stationary-state propagation of Alfvén waves along the thread was derived by
Melis et al. (2021), namely

$\frac{\partial^{2}B_{\perp}}{\partial z^{2}}+\frac{\frac{\partial}{\partial
z}\left(v_{\rm A}^{2}-i\omega\eta_{\rm C}\right)}{\left(v_{\rm
A}^{2}-i\omega\eta_{\rm C}\right)}\frac{\partial B_{\perp}}{\partial
z}+\frac{\omega^{2}}{\left(v_{\rm A}^{2}-i\omega\eta_{\rm
C}\right)}B_{\perp}=0,$ (14)

where $B_{\perp}$ is the transverse perturbation of the magnetic field,
$\omega$ is the angular wave frequency, related with the linear frequency as
$f=\omega/2\pi$, $v_{\rm A}=B/\sqrt{\mu_{0}\rho}$ is the Alfvén speed, with
$\mu_{0}$ the magnetic permeability, and $\eta_{\rm C}$ is the Cowling’s
diffusion coefficient. Again, we note that both $v_{\rm A}$ and $\eta_{\rm C}$
are functions of $z$.

Cowling’s diffusion is the joint effect of Ohm’s and ambipolar diffusion in a
partially ionised plasma. Ohm’s diffusion is caused by the collisions of
electrons with other particles, whereas the ambipolar diffusion is related
with the collisions between neutrals and charged particles. Their coefficients
are written down as

$\displaystyle\eta$ $\displaystyle=$
$\displaystyle\frac{\alpha_{e}}{\mu_{0}e^{2}n_{e}^{2}},$ (15)
$\displaystyle\eta_{\rm A}$ $\displaystyle=$
$\displaystyle\frac{\xi_{n}^{2}}{\mu_{0}\alpha_{n}},$ (16)

where $\eta$ is the Ohm’s coefficient and $\eta_{\rm A}$ is the ambipolar
coefficient, with $\xi_{n}=1-\xi_{i}$ the fraction of neutrals. Then, the
Cowling’s or total diffusion coefficient is computed as $\eta_{\rm
C}=\eta+B^{2}\eta_{\rm A}$.

The wave equation (Equation (14)) needs to be solved considering appropriate
boundary conditions at the ends of the thread. We assume that the wave driver
in located at the left end of the thread, $z=-L/2$. Alfvén waves with a
frequency, $f$, are driven with a prescribed amplitude that depends upon the
adopted spectral weighting function, $A(f)$. Following Tu & Song (2013) and
Arber et al. (2016), the spectral weighting function is assumed to be a piece-
wise power-law function, expressed as

$A\left(f\right)=A_{0}\left\\{\begin{array}[]{lll}\left(\frac{f}{f_{b}}\right)^{5/6},&\textrm{if}&f\leq
f_{b},\\\
\left(\frac{f}{f_{b}}\right)^{-5/6},&\textrm{if}&f>f_{b},\end{array}\right.$
(17)

where $f_{b}$ is the break frequency and $A_{0}$ is a constant. The break
frequency is set to 1.59 mHz as in Tu & Song (2013). The choice of this
frequency is based on the observed spectrum of horizontal photospheric
motions, which suggests that this frequency is between 1 and 10 mHz (see,
e.g., Matsumoto & Kitai, 2010). Physically, this break frequency should
correspond to the beginning of the inertial range governed by the photospheric
turbulence. The constant $A_{0}$ depends on the wave energy flux injected by
the driver. The energy flux for an Alfvén wave of frequency $f$ averaged over
a one full period $1/f$ is

$\langle\@vec{\Pi}\rangle=-\frac{1}{2\mu_{0}}\textrm{Re}\left[v_{\perp}B_{\perp}^{*}\right]\@vec{B},$
(18)

where * denotes the complex conjugate and $v_{\perp}$ is the transverse
velocity perturbation, which can be expressed in terms of $B_{\perp}$ as (see
Melis et al., 2021),

$v_{\perp}=\frac{i}{\omega}\frac{v_{\rm A}^{2}}{B}\frac{\partial
B_{\perp}}{\partial z}.$ (19)

The energy flux can be decomposed as $\langle\@vec{\Pi}\rangle$ =
$\langle\@vec{\Pi}\rangle^{\uparrow}$ \-
$\langle\@vec{\Pi}\rangle^{\downarrow}$, where
$\langle\@vec{\Pi}\rangle^{\uparrow}$ and
$\langle\@vec{\Pi}\rangle^{\downarrow}$ are the parallel and anti-parallel
components of $\langle\@vec{\Pi}\rangle$ with respect to the direction of the
background magnetic field. Their expressions are

$\displaystyle\langle\@vec{\Pi}\rangle^{\uparrow}$ $\displaystyle=$
$\displaystyle\frac{1}{8}\sqrt{\frac{\rho}{\mu_{0}}}Z^{\uparrow}Z^{\uparrow*}\@vec{B},$
(20) $\displaystyle\langle\@vec{\Pi}\rangle^{\downarrow}$ $\displaystyle=$
$\displaystyle\frac{1}{8}\sqrt{\frac{\rho}{\mu_{0}}}Z^{\downarrow}Z^{\downarrow*}\@vec{B},$
(21)

where
$Z^{\uparrow,\downarrow}=v_{\perp}\mp\frac{1}{\sqrt{\mu_{0}\rho}}B_{\perp}$
are the so-called Elsässer variables that represent parallel-propagating and
anti-parallel-propagating Alfvénic disturbances, respectively. At $z=-L/2$,
the quantity $\langle\@vec{\Pi}\rangle^{\uparrow}$ corresponds to the wave
energy flux injected by the driver for one particular frequency, $f$. Thus,
$\Sigma_{f}\langle\boldsymbol{\Pi}\rangle^{\uparrow}$ is the total energy flux
injected by the driver for all the frequencies in the spectrum. The constant
$A_{0}$ is computed by assuming that the total injected energy flux is equal
to a prescribed value, which is hereafter denoted by $\langle\pi\rangle$.

The wave driver is actually located at the photosphere, but the thread model
only includes the coronal part. It is known that the transmission of Alfvén
waves from the photosphere to the corona is heavily influenced by the
conditions in the chromosphere, where reflection and dissipation can equally
be important (see Soler et al., 2017, 2019). To account for the chromospheric
filtering, we use the empirical wave energy transmission coefficient from
Soler et al. (2019), $\mathcal{T}(f,B_{ph})$, which is a function of the wave
frequency, $f$, and the magnetic field strength at the photosphere, $B_{ph}$,
namely

$\begin{split}\mathcal{T}(f,B_{ph})\approx
a_{0}(B_{ph})\frac{1}{\sqrt{2\pi\sigma{(B_{ph})}^{2}}}\exp\left(-\frac{\left(\log_{10}f-\mu(B_{ph})\right)^{2}}{2\sigma{(B_{ph})}^{2}}\right)\\\
\times\left[1+\textrm{erf}\left(\frac{\alpha(B_{ph})}{\sqrt{2}}\frac{\log_{10}f-\mu(B_{ph})}{\sigma(B_{ph})}\right)\right],\end{split}$
(22)

where erf is the error function, and $a_{0}$, $\mu$, $\sigma$ and $\alpha$ are
the amplitude, location, scale and shape parameters, respectively, which are
given in Soler et al. (2019) and depend upon $B_{ph}$. Here, we take
$B_{ph}=100$ G. Therefore, the effective spectral weighting function at the
base of the corona is

$A_{\rm eff.}(f,B_{ph})=A(f)\sqrt{\mathcal{T}(f,B_{ph})}.$ (23)

We note that $A_{\rm eff.}$ is proportional to the square root of
$\mathcal{T}(f,B_{ph})$ because there is quadratic relation between the wave
energy and magnetic field perturbations.

Figure 2 compares the photospheric spectral weighting function, $A(f)$, with
the effective coronal one, $A_{\rm eff.}(f)$. Reflection at the chromosphere
decreases the spectrum amplitude in the low frequency range, while the strong
ion-neutral damping produces the rapid decrease of the amplitude as the
frequency increases. Indeed, the highest frequencies in the spectrum are
effectively suppressed by the chromospheric filtering.

Figure 2: Photospheric spectral weighting function (dash-dotted red line) and
effective coronal spectral weighting function (dashed blue line) as functions
of the wave driver frequency. Note that both axes are in logarithmic scale.
For the purpose of this plot, we set $A_{0}=1$.

Hence, the form of the magnetic field perturbation at the driver location is
prescribed as

$B_{\perp}=A_{\rm
eff.}(f,B_{ph})\exp\left(i\Phi(f)\right)\qquad\textrm{at}\qquad
z=-\frac{L}{2},$ (24)

where $0<\Phi(f)<2\pi$ represents random phase different for each frequency in
the spectrum. On the other hand, the right end of the thread is treated as a
perfectly reflecting boundary for the waves. Contrary to the empirical
transmissivity used at the left boundary, there is no available empirical
reflection coefficient to the used at the right boundary to account for the
effect of the chromosphere, so that a boundary condition has to be imposed
arbitrarily. Melis et al. (2021) tested different boundary conditions, namely
total reflection, total transmission, and partial reflection. Concerning the
wave heating rate in the cool part of the thread, they found no significant
differences between the results for the different boundary conditions. The
physical explanation for this finding is that the high-frequency waves that
actually produce most of the heating have short damping lengths and are almost
completely dissipated during their first passing through the thread, so that
little energy is left to reach and reflect at the boundary. Thus, we restrict
ourselves to the perfectly reflecting scenario, for which the boundary
condition at the right end is

$\frac{\partial B_{\perp}}{\partial z}=0\qquad\textrm{at}\qquad
z=\frac{L}{2}.$ (25)

We discretise the spectrum of the driver into 101 individual frequencies
between 0.1 and 100 mHz with a logarithmic spacing. For each one of the
frequencies in the discretised spectrum, Equation (14) with the above boundary
conditions (Equations (24) and (25)) is numerically solved in _Wolfram
Mathematica_ using again the routine `NDSolve`. However, now we use finite
elements to perform the integration with a customised non-uniform mesh whose
resolution depends upon the considered frequency, $f$. Results from Terradas
et al. (2021) point out that temperature and density profiles are expected to
have sharp variations at the prominence-corona transition region (PCTR).
Concerning the Alfvén waves, the PCTR needs to be resolved with a sufficiently
high resolution to correctly compute the transmission of the waves into the
dense part of the thread. Otherwise, artificial reflections may happen, which
can heavily influence the results. The custom grid is divided into five
different regions: the central part where the plasma is coolest and partially
ionised, two PCTR zones surrounding the centre where there is the sharp
variation of temperature and density, and two outermost parts that represent
the coronal regions. For a prescribed wave frequency, the mesh resolution in
the central and coronal regions is determined by the local Alfvén wavelength
at the centre, $\lambda_{0}=v_{\rm A}(z=0)/f$, or at the ends,
$\lambda_{c}=v_{\rm A}(z=L/2)/f$, of the thread, respectively. We set the
resolution in each region equal to a small fraction of their respective local
wavelengths, being $\lambda_{0}/100$ for the central region and
$\lambda_{c}/20$ for the coronal regions. In the PCTR zones, a fixed mesh
resolution of 100 m is used for all frequencies. The width of the two high-
resolution PCTR zones is set to 40 km and their location is determined by a
routine that detects the position at which the derivative of the temperature
profile with respect to $z$ is maximum or minimum.

The numerical solution of Equation (14) provides us with the magnetic field
perturbation, $B_{\perp}$, along the thread caused by the propagation of
Alfvén waves with frequency, $f$. The next step is to compute the wave heating
rate. The plasma heating is the consequence of the wave energy absorption,
i.e., the part of the wave energy flux that is deposited into the plasma due
to Cowling’s diffusion. The instantaneous heating rate is computed as

$Q_{f}=\mu_{0}\eta_{\rm C}\left|\@vec{j_{\perp}}\right|^{2},$ (26)

where $\@vec{j_{\perp}}$ is the perpendicular component of the current density
given by

$\@vec{j_{\perp}}=\frac{1}{\mu_{0}}\frac{\partial B_{\perp}}{\partial
z}\hat{e}_{\perp},$ (27)

where $\hat{e}_{\perp}$ denotes the unit vector in the direction perpendicular
to the magnetic field. The instantaneous heating rate needs to be averaged
over one full period of the wave, $1/f$, to obtain the net heating produced by
that particular frequency, namely

$\langle Q\rangle_{f}=\frac{\eta_{\rm C}}{2\mu_{0}}\left|\frac{\partial
B_{\perp}}{\partial z}\right|^{2}.$ (28)

In order to obtain the total heating rate, we add together the net heatings
produced by all the frequencies in the broadband spectrum. Thus,

$\langle Q\rangle=\sum_{f}\langle Q\rangle_{f}.$ (29)

The heating profile computed in such a way is the one that is used in the
energy balance equation (Equation (2)).

###  2.3 Self-consistent strategy

The goal is to obtain prominence thread models in which the Alfvén wave
heating obtained from the solution of Equation (14) satisfies the energy
balance condition (Equation (2)). Thus, threads heated by Alfvén waves would
be in thermal equilibrium. This equilibrium is studied in terms of two
parameters: the central temperature, $T_{0}$, and the wave energy flux
injected at the photosphere, $\langle\pi\rangle$. For $T_{0}$, we consider
values between 7,000 to 10,000 K, as consistent with the temperatures in
prominence cores. Regarding the injected energy flux, it is expected that the
wave heating rate increases with the injected flux. Thus, we will
progressively increase the value of $\langle\pi\rangle$ until the requirement
of Equation (13) is no longer satisfied and so no equilibrium models are
possible for such conditions.

The construction of self-consistent thread models heated by Alfvén waves poses
a circular problem. We must solve Equation (14) to compute the wave heating
rate. Since the Alfvén speed and the Cowling’s coefficient appear in Equation
(14), we need the temperature and density profiles to be known beforehand.
However, the temperature profile is computed from Equation (2), which in turn
require the wave heating rate to be known. The solution to this circular
problem is here attacked using an iterative strategy.

To start with, we select particular values of $T_{0}$ and $\langle\pi\rangle$
for which the self-consistent thread model is to be computed. A thread model
satisfying the energy balance condition (Equation (2)) is obtained for the
case with no wave heating, i.e., $\langle Q\rangle=0$. Subsequently, the
Alfvén wave equation (Equation (14)) is solved in that initial thread model
and the wave heating rate is obtained. This heating rate is put back in
Equation (2) and a new thread model, now including wave heating, is computed.
This completes the first iteration of the process. Generally, the initial
thread model obtained for $\langle Q\rangle=0$ and that obtained after the
first iteration would be different. To quantify the difference between models
we study the variation of the length of the cool part of the thread, hereafter
called the ‘thread length’ and denoted by $a$. The thread lengths are computed
by automatically detecting the position of the two PCTR zones, as explained in
Section 2.2. The location of PCTR and so the thread length are computed by
using the second derivative of temperature and detecting the values of $z$ for
which it is zero. We use the parameter $\varepsilon$, which is defined as the
relative error norm of the thread length, namely

$\varepsilon=\frac{|a_{i}-a_{i-1}|}{a_{i}},$ (30)

where $a_{i}$ and $a_{i-1}$ are the thread lengths for the iteration $i$ and
the previous iteration, respectively. We shall assume that a self-consistent
model has been achieved when $\varepsilon<\varepsilon_{0}$, where
$\varepsilon_{0}$ is a prescribed small tolerance. We use
$\varepsilon_{0}=10^{-7}$. Therefore, the process explained before for the
first iteration would be repeated as many times as necessary until converge,
if possible, is reached. If the requirement of Equation (13) is satisfied, the
value of $\varepsilon$ decreases with the number of iterations, leading to
successful convergence. However, if the requirement of Equation (13) is not
satisfied, convergence can never be achieved and the value of $\varepsilon$
may even increase with the number of iterations. Figure 3 displays the
variation of $\varepsilon$ with the number of iterations in two example cases:
a case in which the chosen values of $T_{0}$ and $\langle\pi\rangle$ allow
convergence and a self-consistent model is obtained and another case that does
not converge.

Figure 3: Variation of the parameter $\varepsilon$ with the number of
iterations for a case that converges (dash-dotted red line) and a case that
does not converge (dashed blue line). The tolerance value of
$\varepsilon_{0}=10^{-7}$ is denoted by a horizontal line. Note that the
vertical axis is in logarithmic scale.

Figure 4 summarises in a schematic way how the iterative strategy works to
compute a self-consistent model.

Figure 4: Flux diagram of the iterative strategy that we use to compute self-
consistent thread models. All parameters are defined in the text.

##  3 Results

###  3.1 A typical thread model

We start the presentation of results by analysing a typical thread model
obtained for a certain combination of $T_{0}$ and $\langle\pi\rangle$. This
will help us to understand better the results of the subsequent parameter
study.

Figure 5: Equilibrium profiles along the thread: (a) temperature, (b) density,
(c) ionisation fraction, (d) Ohm’s diffusion coefficient, (e) ambipolar
diffusion coefficient, and (f) volumetric wave heating rate. Results for
$T_{0}=$ 7,000 K and $\langle\pi\rangle=$ 0.2 W m-2.

Hence, we consider $T_{0}=$ 7,000 K and $\langle\pi\rangle=$ 0.2 W m-2. Figure
5 displays the profiles of temperature, density, ionization fraction, Ohm’s
and ambipolar coefficients, and heating rate along the thread. The temperature
profile (Figure 5a) is directly computed from Equation (2) and its shape is
similar to that obtained in the models of Terradas et al. (2021) with a
constant heating term. The temperature profile can be divided into three
distinct regions: the cold central region with a parabolic-like shape around
the temperature minimum, then a sharp increase of the temperature by several
orders of magnitude in a very thin layer that clearly matches the expected
PCTR mentioned before, and, finally, a wider region with a more gentle
increase of temperature that extends up to the ends of the thread, where
typical coronal temperatures of order of $\sim 10^{6}$ K are reached. The
temperature profile is largely symmetric with respect to $z=0$. Very small
asymmetries exist, which are not noticeable at the scale of the figure. In
this particular model, the thread length is $a\approx$ 2.66 Mm, which
corresponds to a small fraction of the total length of the magnetic field
line, namely $a/L\approx 0.027$.

Figure 5b shows the density profile, which is computed from the temperature
profile with the help of Equation (1). The density adopts an inverse profile
with respect to that of the temperature. The largest density of $\sim
10^{-11}$ kg m-3 is located at the centre. In the outermost coronal part of
the model, the density decreases to $\sim 10^{-13}$ kg m-3 near the ends. The
ionisation fraction is plotted in Figure 5c. We recall that the ionisation
fraction is computed from the temperature using Table 1. We obtain a very
narrow partially ionised region around the thread centre, where the plasma is
coldest and densest, with a minimum of $\xi_{i}\approx 0.5$. The plasma
rapidly gets fully ionised in the PCTR, so that $\xi_{i}=1$ in the coronal
part of the model.

The Ohm’s diffusion coefficient, $\eta$, can be seen in Figure 5d.
Essentially, the spatial profile of $\eta$ follows that of the density. The
maximum of $\eta$ is located at the central part, $\eta\sim$
$1\text{\times}{10}^{3}\text{\,}{\mathrm{m}}^{2}\text{\,}{\mathrm{s}}^{-1}$,
and decreases sharply towards the coronal part of the model, where $\eta\sim$
$1\text{\,}{\mathrm{m}}^{2}\text{\,}{\mathrm{s}}^{-1}$. Therefore, Ohm’s
diffusion is three orders of magnitude more efficient in the cool central part
than in the hot coronal part. The ambipolar diffusion coefficient, $\eta_{\rm
A}$, is shown in Figure 5e. The ambipolar coefficient is zero in the coronal
region of the model, where the plasma is fully ionised. Inside the partially
ionised zone, $\eta_{\rm A}$ has two relative maxima slightly displaced from
the centre of the thread. The maximum value is $\eta_{\rm A}\sim$
$1\text{\times}{10}^{13}\text{\,}{\mathrm{m}}^{2}\text{\,}{\mathrm{s}}^{-1}\text{\,}{\mathrm{T}}^{-2}$.
Considering that the background magnetic field strength is $B=10$ G, the
maximum value of the effective ambipolar diffusion coefficient is
$B^{2}\eta_{\rm A}\sim$
$1\text{\times}{10}^{7}\text{\,}{\mathrm{m}}^{2}\text{\,}{\mathrm{s}}^{-1}$,
which points out that ambipolar diffusion is much more efficient than Ohm’s
diffusion in the partially ionised region.

The wave heating rate consistent with this particular model can be seen in
Figure 5f. The heating is very non-uniform along the thread. When approaching
the cold central region, the heating rate increases several orders of
magnitude compared to the values obtained in the hot coronal part. Ambipolar
diffusion is the dominant dissipation mechanism in our model and produces most
of the heating in the cold partially ionised region. In addition, unlike the
rest of quantities, the heating rate displays a clear asymmetry. Towards the
right end, the heating rate tends to zero. The reason for this asymmetry
resides in the fact that the wave driver is only located at the left end of
the model, while the boundary condition at the right end imposes total
reflection. However, the asymmetry is much less noticeable in the central
region and, indeed, the profiles obtained for the rest of model quantities are
nearly symmetric with respect to $z=0$.

###  3.2 Effect of the central temperature

Figure 6: Same as Figure 5, but for different values of the central
temperature, namely $T_{0}=$ 7,000 K (solid blue line), $T_{0}=$ 8,000 K
(dashed red line), $T_{0}=$ 9,000 K (dash-dotted green line), and $T_{0}=$
10,000 K (solid black line). Only a close-up view of the central part and the
PCTR of the thread for $z\geq 0$ is displayed. Results for
$\langle\pi\rangle=$ 0.2 W m-2.

After presenting the results for a particular case, here we explore the effect
of changing the central temperature, $T_{0}$. In all cases, we keep
$\langle\pi\rangle=$ 0.2 W m-2 as before. According to Terradas et al. (2021),
the thread length decreases when the central temperature increases. This
result is confirmed in Figure 6, which displays the equilibrium profiles for
thread models with different values of the central temperature. Due to the
near symmetry of the models, Figure 6 only shows the profiles for positive
values of $z$ corresponding to the central part of the thread and the PCTR,
because it is there where the main differences appear, while the profiles in
the coronal part would be seen as nearly identical.

In agreement with Terradas et al. (2021), when the central temperature
increases, the cold part of the thread becomes narrower. The decrease of the
thread length can be explained by the relative importance of the various terms
in the energy balance equation (Equation (2)). If the effect of heating is
small, which would happen for the small wave energy flux assumed here, thermal
conduction is mostly responsible for balancing radiative losses in the
equilibrium (Terradas et al., 2021). As $T_{0}$ increases, radiative losses
become more and more efficient in the central region. Hence, thermal
conduction requires a stronger temperature gradient to compensate for the
increase of radiation, which results in narrower threads.

Another consequence of the increase of the central temperature is that the
partially ionised plasma gets confined to a thinner region and the lowest
value of $\xi_{i}$ increases. This is explained by the relation between
temperature and ionisation fraction as seen in Table 1. The fact that the
plasma gets more ionised as $T_{0}$ increases impacts on the efficiency of
Ohm’s and ambipolar diffusion. The coefficients associated with both effects
become smaller in magnitude when $T_{0}$ increases.

Despite the fact that the thread length changes, the shapes of the profiles
for the different values of the central temperature are quite similar, except
that for the ambipolar diffusion coefficient (Figure 6e). For high central
temperatures, the maximum of $\eta_{\rm A}$ is located at the centre of the
thread instead of being displaced from the centre when the central temperature
is low. The ambipolar diffusion coefficient has a complicated dependence on
the temperature and ionization degree that results in this peculiar behaviour.
In turn, the heating rate (Figure 6f) has also a similar shape for all the
central temperatures considered, but its value at the thread centre increases
or decreases with $T_{0}$ according to the behaviour of the ambipolar
diffusion coefficient at $z=0$.

###  3.3 Effect of the wave energy flux

Figure 7: Same as Figure 6 with $T_{0}=7000$ K and three different values of
injected wave energy flux, namely $\langle\pi\rangle=$ 0.2 W m-2 (solid blue
line), $\langle\pi\rangle=$ 0.4 W m-2 (dashed red line), and
$\langle\pi\rangle=$ 0.8 W m-2 (dash-dotted green line).

Besides the role of the central temperature, which was already explored by
Terradas et al. (2021), here we are more interested in the effect of the
injected wave energy flux. To this end, we have set the central temperature to
$T_{0}=$ 7,000 K and have considered three different values of the injected
energy flux, namely $\langle\pi\rangle=$ 0.2, 0.4, and 0.8 W m-2. The computed
models are displayed in Figure 7 where, again, only positive values of $z$
corresponding to the central part of the thread and the PCTR are displayed.

Two main effects of varying $\langle\pi\rangle$ are seen in Figure 7. First of
all, the larger $\langle\pi\rangle$, the larger the wave heating rate $\langle
Q\rangle$ (see Figure 7f). This result is not surprising, because the wave
heating rate is expected to be proportional to the injected wave energy flux.
In other words, the larger the available wave energy, the larger the energy
that is dissipated. The second effect that is obvious in Figure 7 is the
growth of the thread length as $\langle\pi\rangle$ increases. Indeed, this
result is a consequence of the increase of the heating rate discussed before.
Following Terradas et al. (2021), we can understand the increase of the thread
length by using an approximate analytic expression for this quantity given by

$a\approx\sqrt{\frac{2\kappa_{\parallel}T_{0}}{\mathcal{L}-\langle
Q\rangle}},$ (31)

where $\kappa_{\parallel}$, $\mathcal{L}$, and $\langle Q\rangle$ need to be
evaluated at the thread centre, $z=0$. If the central temperature remains
constant and the wave heating rate increases at $z=0$ and approaches the
minimum value of radiative loses, then the denominator in Equation (31)
decreases. Consequently, the thread length increases. In Terradas et al.
(2021) a similar conclusion is reached although for a spatially uniform
heating term. Here, the heating rate is very non-uniform along the field line,
but Equation (31) can still be used to understand to effect of the heating on
the thread length. In this regard, and for a fixed central temperature, only
the value of $\langle Q\rangle$ at the centre seems to be relevant, while its
spatial profile does not play a role according to Equation (31). The effect of
the wave heating on the thread length is investigated further in Section 3.5.

###  3.4 Role of wave heating in the energy balance

Figure 8: Comparison of the various terms in the energy balance equation
(Equation (2)) in models with $T_{0}=$ 9,000 K and two values of the injected
wave energy flux: (a) $\langle\pi\rangle=$ 0.1 W m-2 and (b)
$\langle\pi\rangle=$ 30 W m-2. The solid black line corresponds to the wave
heating rate, the dot-dashed blue line denotes radiative losses, and the
dashed red line represents the thermal conduction term. Only a close-up view
of the central part and the PCTR of the thread is displayed.

Here we explore the role of wave heating in the energy balance. We consider an
illustrative example of self-consistent thread models obtained for a central
temperature of 9,000 K and two different injected wave energy fluxes, namely
$\langle\pi\rangle=$ 0.1 and 30 W m-2. Figure 8 compares the importance along
the thread of the various terms that appear in the energy balance equation
(Equation (2)), namely radiation losses, thermal conduction, and wave heating.
To better visualise the results, Figure 8 only displays the cool central part,
the PCTR, and the beginning of the coronal part.

The shape of the wave heating rate is similar for both values of
$\langle\pi\rangle$. The maximum of $\langle Q\rangle$ is located at the
centre of the thread, then there is a sharp decrease of several orders of
magnitudes in the PCTR, so that the values of $\langle Q\rangle$ in the
coronal part of the model are much smaller than in the central cool part.
Indeed, as already pointed out, $\langle Q\rangle$ seems to mimic the profile
of the Cowling diffusion coefficient, $\eta_{\rm C}$, in the central and PCTR
regions, while in the coronal zone it approximately follows the profile of the
Ohm’s diffusion coefficient, $\eta$. This highlights the dominance of each
mechanism in the different regions.

In the case with $\langle\pi\rangle=$ 0.1 W m-2, the wave heating rate is much
less important than radiative losses and thermal conduction, even in the
central cool part. In this case, wave heating has an almost negligible
influence on the energy balance, since radiative losses and thermal conduction
essentially compensate each other everywhere along the field line. Radiative
loses have a minimum located at the centre of the thread, while their maximum
is located inside the PCTR around $T\approx$ 58,000 K. In this case, the
thermal conduction term just follows the spatial dependence of radiative
losses. Comparing the results with the work of Terradas et al. (2021), the
shapes of the radiative losses and thermal conduction terms obtained here are
identical to those obtained in the previous work, although Terradas et al.
(2021) used a constant heating instead of a function depending on position.

The situation changes when we consider a larger wave energy flux of
$\langle\pi\rangle=$ 30 W m-2. Now, a significantly larger heating rate is
obtained in the central region, which approaches the minimum of radiative
losses. For this large value of the wave energy flux, we are close to the
upper boundary of $\langle Q\rangle$ imposed by Equation (13) to obtain a
self-consistent model. Around the thread centre, wave heating compensates a
large fraction of radiation losses, while thermal conduction has now a lower
weight on the energy balance. The relevance of wave heating only applies in
the central part of the thread, where the plasma is cool and partially
ionised, while in the PCTR and, specially, in the coronal region its
contribution remains almost negligible and thermal conduction entirely
compensates cooling. These results are consistent with the work of Melis et
al. (2021), where the heating rate in the cool central region was estimated to
be as important as radiative loses for large enough wave energy fluxes.

As already discussed in Section 3.3, another difference between the results
for the two considered wave energy fluxes is that the length of the thread
increases when $\langle\pi\rangle=$ 30 W m-2 compared to the case with
$\langle\pi\rangle=$ 0.1 W m-2. The computed lengths are $a\approx$ 0.59 Mm
for $\langle\pi\rangle=$ 0.1 W m-2 and $a\approx$ 0.77 Mm for
$\langle\pi\rangle=$ 30 W m-2.

###  3.5 Thread length: parameter survey

In view that the main effect of wave heating is to modify the thread length,
here we perform a parameter study on the variation of this quantity as a
function of $T_{0}$ and $\langle\pi\rangle$. Although some results have
already been obtained in previous sections, now we carry out a more detailed
investigation.

Figure 9: (a) Thread length, $a$, as a function of the injected wave energy
flux, $\langle\pi\rangle$, for different values of central temperature. (b)
Percentage increase of the thread length with respect to the value without
wave heating as a function of $\langle\pi\rangle$ for the same central
temperatures. The vertical lines in panel (b) denote the value of
$\langle\pi\rangle_{\rm max.}$ corresponding to each temperature.

Figure 9a shows the variation of thread length as a function of the injected
wave energy flux for the various central temperatures considered. For each
central temperature, the wave energy flux has been varied between
$\langle\pi\rangle=$ 0.1 W m-2 and the maximum flux that allows the
requirement of Equation (13) to be satisfied, hereafter denoted by
$\langle\pi\rangle_{\rm max.}$. Comparing the results for different
temperatures, we find that for low central temperatures the thread length is
larger than for high central temperatures, as already discussed in Section
3.2, but the evolution of the thread length when $\langle\pi\rangle$ increases
is similar for all the considered central temperatures. Three different
regimes are found. First, for $\langle\pi\rangle\ll\langle\pi\rangle_{\rm
max.}$, the thread length smoothly increases with $\langle\pi\rangle$
following an approximate linear dependence. Then, when $\langle\pi\rangle$
approaches $\langle\pi\rangle_{\rm max.}$, the thread length starts to
increase more abruptly in an approximately exponential fashion. Finally, for
$\langle\pi\rangle\to\langle\pi\rangle_{\rm max.}$, the thread length
asymptotically tends to infinity. Although our iterative strategy allows us to
consider large values of $\langle\pi\rangle$ that differ from
$\langle\pi\rangle_{\rm max.}$ by a tiny percentage, the method becomes
unstable when $\langle\pi\rangle$ is too close to $\langle\pi\rangle_{\rm
max.}$. For this reason, the third regime, i.e., the asymptotic behaviour of
the thread length when $\langle\pi\rangle\to\langle\pi\rangle_{\rm max.}$ is
only partially captured by the results shown in Figure 9a. Such asymptotic
behaviour can be better visualised in Figure 9b, which displays the percentage
increase of the thread length as a function of $\langle\pi\rangle$ with
respect to the length obtained in the absence of wave heating. Depending on
the central temperature, a percentage increase of the thread length between
15% and 30% is obtained for wave energy fluxes just below
$\langle\pi\rangle_{\rm max.}$.

Figure 10: Thread length, $a$, as a function of the central temperature,
$T_{0}$, in the absence of wave heating (lower red line) and for
$\langle\pi\rangle\approx\langle\pi\rangle_{\rm max.}$ (upper blue line). The
grey area indicates the range of possible thread lengths for which a self-
consistent model can be found.

Figure 10 presents the same results of Figure 9 but in a complementary way.
Figure 10 displays the thread length as a function of the central temperature
in two cases: the length obtained when no wave heating is present, i.e., for
$\langle\pi\rangle=0$, and the length obtained for a wave energy flux just
below $\langle\pi\rangle_{\rm max.}$, i.e., at the beginning of the asymptotic
regime. As stated before, the thread length decreases when the central
temperature increases, which is consistent with the previous results of
Terradas et al. (2021). For a given central temperature, Figure 10 shows the
range of values of the thread length for which our iterative method converges
and a self-consistent model is obtained. Threads with smaller lengths are not
possible. Threads with larger lengths are indeed possible, but only in the
asymptotic regime when $\langle\pi\rangle\to\langle\pi\rangle_{\rm max.}$,
which happens in a extremely narrow range of values of $\langle\pi\rangle$
just below $\langle\pi\rangle_{\rm max.}$.

Another interesting result is that $\langle\pi\rangle_{\rm max.}$ is a
function of the central temperature (see Figure 11). High central temperatures
are able to accommodate larger wave energy fluxes than low central
temperatures. Again, the reason for this result can be found in the
requirement of Equation (13) that an equilibrium model must satisfy. At $z=0$,
the radiative losses provided by Athay’s function (Equation 10) increase with
$T_{0}$. Therefore, the maximum wave heating rate at $z=0$ that is consistent
with Equation (13) also becomes larger when $T_{0}$ increases. Obviously,
larger wave heating rates, $\langle Q\rangle$, are associated with larger wave
energy fluxes, $\langle\pi\rangle$. The results of Figure 11 suggest a power-
law dependence of $\langle\pi\rangle_{\rm max.}$ with $T_{0}$. Therefore, we
have performed a least squares fit to the results of Figure 11 considering the
following form,

${\rm log}_{10}\langle\pi\rangle_{\rm max.}=aT_{0}+b,$ (32)

where $a=(770\pm 18)\times 10^{-6}$ and $b=-5.49\pm 0.15$ with
$\langle\pi\rangle_{\rm max.}$ expressed in W m-2 and $T_{0}$ in K. The
$R^{2}$ coefficient of the fit is $R^{2}=0.997$. The fitted function has been
over-plotted in Figure 11, showing a very good approximation to the numerical
data. Hence, for a given central temperature, Equation (32) can be used to
predict the maximum wave energy flux that allows an equilibrium model to
exist.

Figure 11: Maximum wave energy flux injected at the photosphere,
$\langle\pi\rangle_{\rm max.}$, as a function of the thread central
temperature, $T_{0}$. The symbols are the results of the computations, while
the dashed line in the fit of Equation (32). Note that the vertical axis is in
logarithmic scale.

##  4 Conclusions

In this paper we studied 1D models of prominence thin threads that satisfy
energy balance including heating due to Alfvén waves. The present
investigation has been built upon the previous works of Terradas et al. (2021)
and Melis et al. (2021), where equilibrium models and Alfvén wave heating were
studied separately. The computation of the thread models has been done using
an iterative method involving the solution of the energy balance equation
together with the Alfvén wave equation until convergence to a self-consistent
model of a thin thread is achieved. The Alfvén waves are assumed to be
launched from the photosphere by a broadband driver that injects a prescribed
value of the wave energy flux. The waves are then dissipated in the prominence
by Ohm’s and ambipolar diffusion.

The computed thread models are composed of a cold central region that
corresponds to the prominence core (the prominence thread itself), a thin PCTR
in which there is a sharp increase of temperature, and an external region with
solar coronal properties. The wave heating rate is maximum at the central
part, where the plasma is partially ionised and ambipolar diffusion is
responsible for wave dissipation. In the outermost coronal part, wave heating
is negligible because of the irrelevance of Ohm’s diffusion for fully ionised
coronal conditions.

The magnitude of the wave heating rate in the cold central region depends on
the wave energy flux injected by the driver at the photosphere, so that the
larger the injected energy flux, the larger the heating rate. Increasing the
value of the injected energy flux results in obtaining models with longer cold
regions, i.e., longer threads. Self-consistent models can be obtained until
the injected wave energy flux produces a heating rate at the thread centre
that becomes equal to radiative losses. We performed a parametric survey in
order to explore how the thread length correlates with the injected wave
energy flux and the central temperature. The thread length decreases with
increasing central temperature and increases with the injected wave energy
flux. The maximum value of the wave energy flux that allows self-consistent
models to exist has been found to increase with the value of the central
temperature, so that hot threads can accommodate more wave heating than cool
ones.

We have used the injected wave energy flux as a free parameter of the model.
Unfortunately, there are no determinations of the Alfvén wave energy flux
driven at the footpoints of the magnetic structure of prominences. In intense
photospheric flux tubes, calculations of the transverse wave energy flux
driven by horizontal motions show that the flux generated within the
photospheric flux tubes can be as large as $\sim 10^{6}$ W m-2 (see, e.g.,
Spruit, 1981; Ulmschneider, 2000; Noble et al., 2003; Shelyag et al., 2011).
When averaged over the entire photosphere and considering the photospheric
filling factor, the resulting averaged driven flux is of the order of $\sim
10^{4}$ W m-2 (see, e.g., De Pontieu et al., 2001; Goodman, 2011; Tu & Song,
2013; Arber et al., 2016). This estimated wave energy flux applies in the case
of photospheric bright points, but it is probably much larger than the wave
energy flux driven in the case of prominences. If this large energy flux were
applicable to prominences, then our results indicate that no equilibrium
models of threads would be possible: wave heating would be simply too large.
Stable threads could not sustain such large wave energy fluxes. In essence,
the existence or absence of equilibrium models depends upon the balance
between wave heating and radiation losses at the thread centre. While
radiation losses mainly depend on the local properties of the thread, the
magnitude of wave heating is determined by how the waves are driven at the
photospheric regions where the prominence magnetic field is anchored.
Recently, Li et al. (2022) reported Alfvénic waves along prominence threads
with an estimated energy flux of 16.2–37.3 W m-2, which is more compatible
with the wave energy fluxes needed to obtain equilibrium models.

For simplicity and to consider the same set-up as in Melis et al. (2021), a
perfectly reflecting boundary has been implemented at the right end of the
thread instead of considering the presence of two different drivers located at
both ends. The expected results of adding a second driver to the model are
straightforward. If the second driver injects the same energy flux as the
first one, the heating rate in the cool region would double its value. In
turn, the maximum value of the wave energy flux to obtain a self-consistent
model would then need to account for the sum of the fluxes injected by both
drivers.

A limitation of this work is that the computation of the wave heating rate
assumes the stationary state of wave propagation. In other words, it has been
assumed that the waves have been propagating and reflecting along the thread
for enough time to reach a stationary pattern. An alternative approach would
be to solve the time-dependent problem, for which we should resort to full
numerical simulations. The comparison of time-dependent results with those of
the stationary problem is relevant and may be tackled in the future.

Although the present method of constructing thread models heated by Alfvén
waves provides interesting results, the background configuration is still
simple compared to the observed complexity of prominences. A further extension
would be to consider 2D models in cylindrical geometry, which should include
the perpendicular component of thermal conduction. Concerning the Alfvén
waves, the use of 2D models would introduce the presence of resonances in the
Alfvén continuum. Adding these new ingredients would allow us to obtain more
accurate information about the efficiency of wave heating in solar prominences
and the conditions for the existence of equilibrium models. This could also be
undertaken in a follow-up work.

Finally, we connect with the comment made in the Introduction about the
possible role that Alfvén wave heating may have to complement radiative
heating in prominences. Here we have computed thread models with realistic
lengths and central temperatures without including radiative heating. This
suggests that Alfvén wave heating may play not only a complementary role but
an important part in the energy balance of prominence threads that future
works should keep investigating. Ideally, future models should include both
radiative and wave heating, which is a challenging task.

###### Acknowledgements.

This publication is part of the R+D+i project PID2020- 112791GB-I00, financed
by MCIN/AEI/10.13039/501100011033. This research was supported by the
International Space Science Institute (ISSI) in Bern, through ISSI
International Team project #457 (The Role of Partial Ionization in the
Formation, Dynamics and Stability of Solar Prominences). We thank Ramón Oliver
for some comments regarding the importance of correctly resolving the PCTR. We
acknowledge the unknown referee for very constructive comments that helped us
to improve the paper.

## References

  * Anzer & Heinzel (1999) Anzer, U. & Heinzel, P. 1999, A&A, 349, 974
  * Arber et al. (2016) Arber, T. D., Brady, C. S., & Shelyag, S. 2016, ApJ, 817, 94
  * Arregui et al. (2018) Arregui, I., Oliver, R., & Ballester, J. L. 2018, Living Reviews in Solar Physics, 15, 3
  * Athay (1986) Athay, R. G. 1986, ApJ, 308, 975
  * Ballester (2015) Ballester, J. L. 2015, in Astrophysics and Space Science Library, Vol. 415, Solar Prominences, ed. J.-C. Vial & O. Engvold, 259
  * Ballester et al. (2018) Ballester, J. L., Alexeev, I., Collados, M., et al. 2018, Space Sci. Rev., 214, 58
  * Brughmans et al. (2022) Brughmans, N., Jenkins, J. M., & Keppens, R. 2022, A&A, 668, A47
  * Dalgarno & McCray (1972) Dalgarno, A. & McCray, R. A. 1972, ARA&A, 10, 375
  * De Pontieu et al. (2001) De Pontieu, B., Martens, P. C. H., & Hudson, H. S. 2001, ApJ, 558, 859
  * Degenhardt & Deinzer (1993) Degenhardt, U. & Deinzer, W. 1993, A&A, 278, 288
  * Gilbert (2015) Gilbert, H. 2015, in Astrophysics and Space Science Library, Vol. 415, Solar Prominences, ed. J.-C. Vial & O. Engvold, 157
  * Goodman (2011) Goodman, M. L. 2011, ApJ, 735, 45
  * Heasley & Mihalas (1976) Heasley, J. N. & Mihalas, D. 1976, ApJ, 205, 273
  * Heinzel (2015) Heinzel, P. 2015, in Astrophysics and Space Science Library, Vol. 415, Solar Prominences, ed. J.-C. Vial & O. Engvold, 103
  * Heinzel & Anzer (2012) Heinzel, P. & Anzer, U. 2012, A&A, 539, A49
  * Heinzel et al. (2010) Heinzel, P., Anzer, U., & Gunár, S. 2010, Mem. Soc. Astron. Italiana, 81, 654
  * Heinzel et al. (2015) Heinzel, P., Gunár, S., & Anzer, U. 2015, A&A, 579, A16
  * Hermans & Keppens (2021) Hermans, J. & Keppens, R. 2021, A&A, 655, A36
  * Hillier et al. (2013) Hillier, A., Morton, R. J., & Erdélyi, R. 2013, ApJ, 779, L16
  * Li et al. (2022) Li, D., Xue, J., Yuan, D., & Ning, Z. 2022, Science China Physics, Mechanics, and Astronomy, 65, 239611
  * Lin (2011) Lin, Y. 2011, Space Sci. Rev., 158, 237
  * Lin et al. (2009) Lin, Y., Soler, R., Engvold, O., et al. 2009, ApJ, 704, 870
  * Martin (2015) Martin, S. F. 2015, in Astrophysics and Space Science Library, Vol. 415, Solar Prominences, ed. J.-C. Vial & O. Engvold, 205
  * Matsumoto & Kitai (2010) Matsumoto, T. & Kitai, R. 2010, ApJ, 716, L19
  * Melis et al. (2021) Melis, L., Soler, R., & Ballester, J. L. 2021, A&A, 650, A45
  * Mok et al. (1990) Mok, Y., Drake, J. F., Schnack, D. D., & van Hoven, G. 1990, ApJ, 359, 228
  * Noble et al. (2003) Noble, M. W., Musielak, Z. E., & Ulmschneider, P. 2003, A&A, 409, 1085
  * Ofman et al. (1998) Ofman, L., Klimchuk, J. A., & Davila, J. M. 1998, ApJ, 493, 474
  * Parenti (2014) Parenti, S. 2014, Living Reviews in Solar Physics, 11, 1
  * Parenti & Vial (2007) Parenti, S. & Vial, J. C. 2007, A&A, 469, 1109
  * Pécseli & Engvold (2000) Pécseli, H. & Engvold, O. 2000, Sol. Phys., 194, 73
  * Schure et al. (2009) Schure, K. M., Kosenko, D., Kaastra, J. S., Keppens, R., & Vink, J. 2009, A&A, 508, 751
  * Shelyag et al. (2011) Shelyag, S., Fedun, V., Keenan, F. P., Erdélyi, R., & Mathioudakis, M. 2011, Annales Geophysicae, 29, 883
  * Soler et al. (2016) Soler, R., Terradas, J., Oliver, R., & Ballester, J. L. 2016, A&A, 592, A28
  * Soler et al. (2017) Soler, R., Terradas, J., Oliver, R., & Ballester, J. L. 2017, ApJ, 840, 20
  * Soler et al. (2019) Soler, R., Terradas, J., Oliver, R., & Ballester, J. L. 2019, ApJ, 871, 3
  * Soler et al. (2021) Soler, R., Terradas, J., Oliver, R., & Ballester, J. L. 2021, ApJ, 909, 190
  * Spruit (1981) Spruit, H. C. 1981, A&A, 98, 155
  * Terradas et al. (2021) Terradas, J., Luna, M., Soler, R., et al. 2021, A&A, 653, A95
  * Tu & Song (2013) Tu, J. & Song, P. 2013, ApJ, 777, 53
  * Ulmschneider (2000) Ulmschneider, P. 2000, Chromosphere: Heating Mechanisms, Encyclopedia of Astronomy and Astrophysics, ed. P. Murdin (Nature Publishing Group), 2260
  * Vial & Engvold (2015) Vial, J.-C. & Engvold, O. 2015, Solar Prominences, Vol. 415


hlvtwo nohl end inline

# Patterning of ultrathin YBCO nanowires using a new focused-ion-beam process

N Curtz1,2, E Koller2, H Zbinden1, M Decroux2, L Antognazza2, Ø Fischer2 and N
Gisin1 1 Group of Applied Physics, University of Geneva, 1211, Geneva 4,
Switzerland  2 Department of Condensed Matter Physics, University of Geneva,
1211, Geneva 4, Switzerland<EMAIL_ADDRESS>
###### Abstract

Manufacturing superconducting circuits out of ultrathin films is a challenging
task when it comes to pattern complex compounds, which are likely to be
deteriorated by the patterning process. With the purpose of developing high-Tc
superconducting photon detectors, we designed a novel route to pattern
ultrathin YBCO films down to the nanometric scale. We believe that our method,
based on a specific use of a focused ion beam, consists in locally implanting
Ga${}^{\text{3+}}$ ions and/or defects instead of etching the film. This
protocol could be of interest to engineer high-Tc superconducting devices
(SQUIDS, SIS/SIN junctions, Josephson junctions), as well as to treat other
sensitive compounds.

###### pacs:

85.25.Am, 81.16.Nd, 74.78.Bz, 81.15.Cd, 85.25.Pb

##  1 Introduction

Superconducting Single-Photon Detectors (SSPDs) are superconducting devices
designed with the purpose of detecting light down to single-photon level. They
present a good quantum efficiency (QE > 10%), low dark-count rate (DK < 10
Hz), and high operating frequency (> GHz), outperforming in a number of cases
InGaAs Avalanche PhotoDiodes (APDs) [1]. These characteristics make them a
premium candidate for single-photon telecommunication and applications like
Quantum Key Distribution. Their underlying mechanism is based upon the
formation of a hotspot in a current-biased superconducting stripe [2]. The
creation of the hotspot is triggered by the incoming photon whose energy
locally thermalizes the stripe, confining the bias current hence raising its
density, up to the point of overcoming its critical value, resulting in a
local transition of the stripe.

To achieve that, the circuit geometry has to fulfill drastic geometrical
constraints. First, the stripe has to be narrow enough, otherwise the
variation of the current density isn’t sufficiently significant, preventing
the transition to take place and the voltage pulse to be detectable. The
section of the nanowire should also be extremely homogeneous, since any
constriction locally lowers the critical current, hence affects the whole
device with for aftermath a drop of the QE; the closer to 1 is the ratio
$I_{\text{bias}}/I_{\text{c}}$, the closer to the dissipative state is the
device, and therefore the more important is the section homogeneity. Finally,
the detector’s recovery time is governed by the device thickness,
necessitating devices less than 15 nm thick. Whereas such devices have been
successfully produced with low-Tc superconductors such as NbN operated at 2.6
K [3], their realization with high-Tc compounds remains a challenge. High-Tc
SSPDs would nevertheless allow a higher working temperature, hence a
significant reduction of the associated cryogenic costs.

Among high-T${}_{\text{c}}$ materials, cuprates present the advantage of a low
kinetic inductance, leading to fast response times. From a purely structural
point of view,
Nd${}_{\text{1+}x}$Ba${}_{\text{2-}x}$Cu${}_{\text{3}}$O${}_{\text{7+}\delta}$
presents excellent crystallographic and planeity properties [4], which are
interesting features given the aforementioned geometrical constraints of
SSPDs. It turns out however that the intrinsic loose stoichiometry of Nd
atoms, who interdiffuse with Ba ones, leads to a high unstability of the
oxygen content and therefore to an important loss of T${}_{\text{c}}$ during
the patterning process.
YBa${}_{\text{2}}$Cu${}_{\text{3}}$O${}_{\text{7-}\delta}$ (YBCO) is much more
stable and appears to be a good candidate for high-T${}_{\text{c}}$ detectors.
Several routes have been carried out to create junctions or patterned
structures out of high-T${}_{\text{c}}$ thin films, such as Selective
Epitaxial Growth [5], Electron Beam Lithography / Ion Beam Etching (EBL/IBE)
[6], Ion Irradiation [7], using an Atomic Force Microscope [8] or a Focused
Electron Beam Irradiation [9]. Nanobridges have been fabricated with a Focused
Ion Beam [10]. However those experiments were performed on films with
thickness $d$ > 20 nm. Here we report a new method using such an apparatus to
write an arbitrary pattern upon an ultrathin (< 20 nm) film, allowing to
manufacture YBCO superconducting circuits. The key point of this method is
that to produce the structure the superconducting phase is locally altered
rather than etched.

##  2 Experimental

###  2.1 Overview of the FIB-based protocol

We designed a 2-step protocol, involving a preliminary chemical etching
followed by a focused ion beam (FIB) managed nanostructuration. An overview of
the scheme used to create YBCO circuits embedding paths to characterize their
transport properties with 4-point measurements is given in figure 1.

Figure 1:  Overview of the patterning protocol. (a) Photograph of the
preliminary structure created from a thin film by photolithography and
chemical etching. (b) SEM micrograph of a 20 $\mu$m stripe engineered with a
high-current Focused Ion Beam. (c) SEM micrograph of a 12 nm thick, 15 $\mu$m
long, 500 nm wide stripe patterned with a lower current. (d) SEM micrograph of
a 1 $\mu$m wide meandering circuit written with a lower current.

###  2.2 Growth of YBCO films

$c$-axis YBCO films were deposited by RF magnetron sputtering on (100)
SrTiO${}_{\text{3}}$ substrates. Samples are heated to 700°C and exposed to
the sputtering of a YBCO target in a Argon plasma. Along with the Ar flow a
complementary O${}_{\text{2}}$ flow (1:5 ratio) accounts for the growth of a
tetragonal, non-superconducting phase, under a total pressure
P=8.10${}^{\text{-2}}$ mbar. An in-situ 2-hour-long annealing at 580°C in a
600 mbar O${}_{\text{2}}$ atmosphere causes this
YBa${}_{\text{2}}$Cu${}_{\text{3}}$O6 phase to undergo a tetragonal-
orthorhombic transition to optimally-doped superconducting
YBa${}_{\text{2}}$Cu${}_{\text{3}}$O${}_{\text{7-}\delta}$. The critical
temperature of bulk YBCO is 92 K; this critical temperature decreases with the
film’s thickness, down to T${}_{\text{c0}}$(d=12 nm) $\approx$ 80 K. The high
crystallinity of the films was demonstrated with X-ray diffraction
measurements such as depicted in figure 2. Up to 4 degrees the grazing
incidence scan leads to Kiessig fringes contributions from both films [11].
Around the (001) YBCO Bragg peak are located secondary fringes clearly showing
a finite size effect (Laue oscillations), demonstrating the high quality of
the crystallographic layers. The fringes of X-ray spectra allow to determine a
film’s thickness with unit cell accuracy. In the following, all the films
processed are at d=12 nm, with T${}_{\text{c0}}\approx$ 80 K.

Figure 2:  Typical X-ray $\theta$-2$\theta$-spectrum of a 12 nm thick YBCO
film passivated with a 8 nm thick amorphous PBCO layer.

The sputtering process also deposits CuO${}_{\text{2}}$ surface particles
located on the top of the YBCO phase, with a typical diameter estimated at 250
nm; those might indeed be the origin of shortcuts after patterning. No
reproducible way of getting particle-free samples or eliminating them has been
found, however most of the time these particles are not a real obstacle to the
patterning process, since the core structure of the devices can most of the
times be chosen to be located in a clear area (figure 1c). Even when it’s not
the case, such as in figure 1d, table 2 shows that ultimately these particles
are not an issue.

In order to improve electrical contacts, we adapted a gold evaporation system
involving a mechanical mask to deposit gold slots in-situ immediatly after the
YBCO deposition. In addition it was observed that the in-situ deposition of a
8 nm thick amorphous PrBaCuO passivation cap layer over the whole sample
attenuates the loss of critical temperature occurring during the chemical
etching process [12]. Therefore, we end up with a YBCO/Au/PBCO topology.

###  2.3 Photolithography and chemical etching

The patterning of a 50 $\mu$m wide stripe is done by optical lithography and
chemical etching with orthophosphoric acid H${}_{\text{3}}$PO${}_{\text{4}}$
(1%). Four independent structures are patterned in a single run (figure 1a).
One of them, top-left on the figure, is specifically designed to electrically
characterize the sample at this point of the protocol.

###  2.4 FIB writing

As the second step of the protocol, a dual beam (FEI Nova 600 Nanolab scanning
electron microscope / 30 kV focused Ga${}^{\text{3+}}$ ion beam) was used to
carry out a finer patterning of the YBCO films into meandering stripes
suitable for optical characterization of SSPDs. In the context of YBCO thin
layers, some precautions to make a relevant use of the FIB are necessary since
the standard manner of working is too destructive for the samples.
Accordingly, we devised a specific modus operandi to satisfy the needed
requirements.

First of all, due to the extreme thickness of the involved films, the
alignment routines cannot be handled with the standard procedure, as the whole
window, containing areas destined to remain superconducting is exposed to the
Ga${}^{\text{3+}}$ beam during the operation, and irreversibly damaged. To
circumvent this problem we have to synchronize the SEM and the FIB on a non-
critical area, then align the pattern to etch with the sample using the SEM.

Secondly, to find which parameters (intensity, number of passes, numerical
aperture…) are optimal to produce an effective pattern, we tried to separate
two YBCO areas by creating FIB-made barriers and injecting current through
them; results for a 50 pA beam are shown in table 1. We determined that at 50
pA, with a beam diameter $\delta\approx$ 20 nm, 50 passes are optimal to
electrically isolate both areas, while at 1 nA, with a beam diameter
$\delta\approx$ 100 nm, 10 passes are optimal. The fact that for N>20 the
resistance exceeds 3 M$\Omega$ rules out an exclusive YBCO${}_{\text{7}}$ to
YBCO${}_{\text{6}}$ transformation due to oxygen loss, for if it were the
case, the barrier transition would have resulted in a 5 k$\Omega$ resistance
(assuming very conservatively $\rho$(underdoped)=10${}^{\text{-1}}$
$\Omega$.cm at room temperature [13]) for a 1 $\mu$m wide affected width. Even
if the whole stripe had been transformed to YBCO${}_{\text{6}}$ its resistance
wouldn’t have exceeded 150 k$\Omega$. We attribute this result to the fact
that some Ga${}^{\text{3+}}$ ions are implanted or give rise to columnar
defects into the YBCO phase during the exposure, turning it locally into an
insulator.

Table 1: Room-temperature resistance of a FIB-made barrier across a 20 $\mu$m
wide stripe for various numbers of passes (N) with a 50 pA beam.

  * N | 0 | 1 | 2 | 5 | 20 | 50 | 100
---|---|---|---|---|---|---|---
R ($\Omega$) | 2k | 2k | 3k | 700k | 3M | >10M | >10M

Lateral contamination is a crucial parameter to consider since it directly
addresses both the issues of smallest reachable dimensions and current
homogeneity. To illustrate this point, we obtained different kinds of
superconducting transitions on structures of different widths by varying the
number of FIB passes N. Figure 3 shows a comparison between resistive profiles
for 20 $\mu$m wide stripes patterned with respectively N=500 and N=10 passes.
A significant loss of T${}_{\text{c}}$ occurs for N=500. This T${}_{\text{c}}$
loss may be due to a reduction of the O${}_{\text{2}}$ stoichiometry, but in
this case it can’t be attributed to local heating occurring during the writing
process. Indeed, the elevation of the film temperature has been estimated
within a simple thermal model and found to be lower than 5 K, under the beam
writing conditions (voltage: 30 kV, current: 1 nA, diameter: 100 nm). Besides,
it is unlikely that Ga ions can travel 10 $\mu$m accross the sample resulting
in implantations responsible for the T${}_{\text{c}}$ losses observed in
figure 3. Consequently N=500 was ruled out for our process. For N=10 no such
loss happens, although the current is confined in the stripe.

Figure 3:  Resistive curves of 12 nm thick micron-scale circuits patterned
using a Focused Ion Beam with a beam current i=1 nA. N is the number of passes
used to produce the structure.

Finally, to pattern a meandering wire, we first produce a 20 $\mu$m wide
stripe (figure 1b), then we force the current to follow a serpentine path by
adding some crossing barriers (figure 1d). The supporting stripe has to be
manufactured at 1nA because of window size considerations. Thus we choose to
use two different FIB currents to pattern a meander : 1 nA, 10 passes for the
supporting stripe; 50 pA, 50 passes for the meander itself. Supplementary
barriers parallel to the stripe are added to avoid side effects near to the
1nA-exposed areas. Using this protocol superconducting straight stripes down
to 500 nm wide (figure 1c), and meanders down to 1 $\mu$m wide (figure 1d)
were produced. It’s worth mentioning that at no point we were confronted to
some limitation so in principle it should be possible to achieve even better
resolutions with this method.

##  3 Results and discussion

The $\rho$ vs T curves of the samples can be followed during the different
steps of the protocol: after the film deposition the resistivity is measured
with the Van der Pauw method [14]; after chemical etching a four-point
measurement is carried out along a 50 $\mu$m wide stripe using the top-left
structure of figure 1a. Eventually the final circuit is characterized as shown
in figure 1b. Results are plotted in figure 4, and demonstrate that the
obtained circuits are superconducting, although a small loss of
T${}_{\text{c}}$ as well as a broadening of the resistive transition is
observed. This could be explained by the fact that the samples are
intrinsically inhomogeneous, and by restricting the superconducting geometry
to narrow circuits new areas with lower T${}_{\text{c}}$ enter the current
path; but it could also be material damages occuring during the processing. On
another hand, $\rho_{\text{100K}}$ and $\rho_{\text{300K}}$ are slightly
higher at the end of the process.

As mentioned previously the writing routine with the focused ion beam doesn’t
etch the film; instead the beam turns locally the irradiated area into an
insulating phase by implanting columnar defaults inside it. The penetration
depth of Ga${}^{\text{3+}}$ ions is about 70 nm, far superior to the YBCO and
amorphous PBCO combined thickness. From this point of view this method differs
from EBL/IBE processes where parts of the film are physically removed to
create the superconducting pattern, but it also implies caution if one desires
to implement it on with thicker films. This approach prevents the YBCO phase
from being in contact with air, which could be an escape path for oxygen, both
with the PBCO top passivation layer and the lateral insulating phase. Table 2
presents resistance measurements obtained for different structures, the last
column showing the consistency of the results and ensuring that the method
presented in the paper leads to reproducible structures.

Figure 4:  Resistivity vs temperature for different samples, with additional
curves at intermediary steps of the process.  Table 2: Summary of
T${}_{\text{c}}$ and resistance for different kinds of samples.
R${}_{\text{norm}}$ is the equivalent resistance computed with the structure
normalized to a 1 $\mu$m wide, 15 $\mu$m long line.

  * shape | length | width | T${}_{\text{c0}}$ | R(100K) | R${}_{\text{norm}}$ ($k\Omega$)
---|---|---|---|---|---
meander | 180 $\mu$m | 2 $\mu$m | 72 K | 32 k$\Omega$ | 5.3
meander | 430 $\mu$m | 1 $\mu$m | 70 K | 190 k$\Omega$ | 6.6
meander | 430 $\mu$m | 1 $\mu$m | 75 K | 140 k$\Omega$ | 4.9
stripe | 15 $\mu$m | 1 $\mu$m | 70 K | 5 k$\Omega$ | 5
stripe | 15 $\mu$m | .8 $\mu$m | 75 K | 6.6 k$\Omega$ | 5.3
stripe | 15 $\mu$m | .5 $\mu$m | 65 K | 12 k$\Omega$ | 6

Figure 5 present the current-voltage and the $\rho$-$j$ characteristics of one
sample. Both present a negative curvature in the whole range of temperature
covered, yielding the existence, in every case, of a true critical current
density $j_{\text{c}}$. This ensures that the superconducting phase isn’t in a
flux-creep state in spite of the thinness of the sample, which is consistent
with the observation that above $d$ = 10 nm flux line lattices behave like a
3D system [15].

The critical current is generally defined using the standard criterion of an
electric field of 1$\mu$V/cm. In our case the measurement noise due to the
high impedance of the line sets a resolution limit of 2 $\mu$V across the
whole meander, corresponding to an equivalent electric field E = 110 $\mu$V/cm
as shown on figures 5a and 5b. We clearly see on figure 5b, especially at T =
65 K, that at 110 $\mu$V/cm (or 2 $\mu$V) the critical current densities are
overestimated and that such a determination isn’t relevant.

However, it has been reported [16] that a precise measurement of the flux flow
resistivity at high current density allows to extract $j_{\text{c}}$ from the
best power law fit $V\propto(j-j_{\text{c}})^{n}$, $j_{\text{c}}$ and $n$
being the fitting parameters. Figure 6 shows the very good fits obtained with
power laws, giving good confidence for the determination of $j_{\text{c}}$’s
by this procedure. Moreover, $n$ is found to be around 5 with little variation
over the whole temperature range. A significant advantage of this fitting
method is that the knowledge of the current-voltage characteristic at high
current density is sufficient to infer $j_{\text{c}}$.

The temperature-dependence of $j_{\text{c}}$ is described within the Ginzburg-
Landau theory of superconductivity as $j_{\text{c}}$ =
$j_{\text{c0}}(1-(\frac{T}{T_{\text{c}}})^{2})(1-(\frac{T}{T_{\text{c}}})^{4})^{1/2}$
[17]. Figure 7 shows that this model fits extremely well our experimental data
and allows to extrapolate $j_{c}(0K)\approx$ 4.1 MA/cm${}^{\text{2}}$, which
is two orders of magnitude smaller than the depairing limit
$j_{\text{d}}=0.54\frac{B_{\text{c}}}{\mu_{0}\lambda}\approx 3~{}10^{8}$
A/cm${}^{\text{2}}$, demonstrating the good quality of the samples after the
processing. It’s worth noting that the hypothesis of inhomogeneities in the
sample previously mentioned is reinforced by the fact that the best Ginzburg-
Landau fit is obtained for T${}_{\text{c}}$ = 68 K, below the critical
temperature of 72 K found with resistivity measurements (see table 2).

A last point concerns long-term behavior of the nanowires. Although no
systematic study was carried out, some samples were characterized several
times at low temperature, with more than one month gone by between the
measurements. No deterioration of any kind was observed, which is an
indication that the manufactured structures are unaffected by time and thermal
cycling.

Figure 5:  (a) Voltage vs current density at various temperatures for a 2
$\mu$m wide meandering circuit. (b) Resistivity vs current density for the
same sample.

Figure 6:  Same data as in figure 5a represented as voltage vs
$j/j_{\text{c}}-1$. Power law fits $V=B(j-j_{\text{c}})^{n}$ (continuous
lines) correctly describe the curves and allow to determine $j_{\text{c}}$. In
inset is reported $n$ as temperature drops from 65 K to 30 K.  Figure 7:
Critical current density vs reduced temperature $t=T/T_{\text{c}}$ for the
same meandering circuit. $j_{\text{c}}$ is obtained through the power law fits
presented in figure 6. The solid line is the best fit with the theoretical
Ginzburg-Landau model.

##  4 Conclusion

To summarize, a modus operandi to create ultrathin superconducting YBCO
circuits by implanting Ga${}^{\text{3+}}$ ions with a Focused Ion Beam was
devised. We confined current in straight stripes down to 500 nm wide and
produced 1 $\mu$m wide meandering wires using a 50 pA beam current. The
consistency of the resistivity vs temperature profiles measured on the samples
at the different steps of the processing ensures the reproducibility of this
patterning method for superconducting films. For one sample, the critical
current density extrapolated to 0 K has been found to be only two orders of
magnitude times smaller than the depairing limit, demonstrating its quality.
The most natural application of this protocol would be the manufacturing of
high-T${}_{\text{c}}$ superconducting devices. Photoresponse experiments to
characterize the devices as single-photon detectors are left for future work.

The authors would like to thank Michaël Pavius, Kevin Lister, Samuel Clabecq,
and Philippe Flückiger for providing access and training to EPFL’s focused ion
beam, as well as Jean-Claude Villégier for helpful discussions. This work is
supported by the European project Sinphonia (contract No. NMP4-CT-2005-16433),
and the Swiss poles NCCR MaNEP and NCCR Quantum Photonics.

## References

  * [1] R T Thew, N Curtz, P Eraerds, W Walenta, J-D Gautier, E Koller, J Zhang, N Gisin, and H Zbinden.  Approaches to single photon detection.  Nucl. Instr. and Meth. A, 610 (1):16–19, 2009.  doi:10.1016/j.nima.2009.05.031.
  * [2] G Gol’tsman, O Okunev, G Chulkova, A Lipatov, A Semenov, K Smirnov, B Voronov, A Dzardanov, C Williams, and R Sobolewski.  Picosecond superconducting single-photon optical detector.  Appl. Phys. Lett., 79:705, 2001.
  * [3] F Marsili, D Bitauld, A Gaggero, R Leoni, F Mattioli, S Hold, M Benkahoul, F Lévy, and A Fiore.  High efficiency NbN nanowire superconducting single photon detectors fabricated on MgO substrates from a low temperature process.  Optics Express, 16:3191, 2008.
  * [4] M Badaye, J G Wen, K Fukushima, N Koshizuka, T Morishita, T Nishimura, and Y Kido.  Superior properties of NdBa${}_{\text{2}}$Cu${}_{\text{3}}$Oy over YBa${}_{\text{2}}$Cu${}_{\text{3}}$Ox thin films.  Supercond. Sci. Technol., 10:825–830, 1997.
  * [5] C A J Damen, H-J H Smilde, D H A Blank, and H Rogalla.  Selective epitaxial growth for YBCO thin films.  Supercond. Sci. Technol., 11:437–448, 1998.
  * [6] H Schneidewind, F Schmidl, S Linzen, and P Seidel.  The possibilities and limitations of ion-beam etching of YBa${}_{\text{2}}$Cu${}_{\text{3}}$O${}_{\text{7-}x}$ thin films and microbridges.  Physica C, 250:191–201, 1995.
  * [7] N Bergeal, X Grison, J Lesueur, G Faini, M Aprili, and J-P Contour.  High Tc superconducting quantum interferences devices made by ion irradiation.  Appl. Phys. Lett., 89:112515, 2006.
  * [8] C Delacour, J Claudon, J-Ph Poizat, B Pannetier, V Bouchiat, R Espiau de Lamaëstre, J-C Villégier, M Tarkhov, A Korneev, B Voronov, and G Gol’tsman.  Superconducting single photon detectors made by local oxydation with an atomic force microscope.  Appl. Phys. Lett., 90:191116, 2007.
  * [9] W E Booij, A J Pauza, E J Tarte, D F Moore, and M G Blamire.  Proximity coupling in high-Tc Josephson junctions produced by focused electron beam irradiation.  Phys. Rev. B, 455:14600, 1997.
  * [10] S-G Lee, S Oh, C S Kang, and S-J Kim.  Superconducting nanobridge made from YBa${}_{\text{2}}$Cu${}_{\text{3}}$O${}_{\text{7}}$ film by using focused ion beam.  Physica C, 460:1468–1469, 2007.
  * [11] K N Stoev and K Sakurai.  Review on grazing incidence X-ray spectrometry and reflectometry.  Spectrochimica Acta Part B, 54:41–82, 1999.
  * [12] A Jaeger, J-C Villégier, P Bernstein, J Bok, and L Force.  Field effect in multilayered YBaCuO/PrBaCuO devices.  Ieee Trans. Appl. Supercond., 3:2933, 1993.
  * [13] K Semba, A Matsudaa, and M Mukaida.  Carrier-concentration-driven superconductor-to-insulator transition in YBa${}_{\text{2}}$Cu${}_{\text{3}}$O${}_{\text{6+}x}$.  Physica B, 281:904–905, 2000.
  * [14] L J Van der Pauw.  A method of measuring the resistivity and Hall coefficient on lamellae of arbitrary shape.  Philips Technical Review, 20:220–224, 1958.
  * [15] J-M Triscone and Ø Fischer.  Superlattices of high-temperature superconductors.  Rep. Prog. Phys., 60:1673–1721, 1997.
  * [16] L Antognazza, M Decroux, S Reymond, E de Chambrier, J-M Triscone, W Paul, M Chen, and Ø Fischer.  Simulation of the behavior of superconducting YBCO lines at high current densities.  Physica C, 372-376:1684–1687, 2002.
  * [17] C P Poole, H A Farach, and R J Creswick.  Superconductivity.  Academic Press, 1995.