# Discovery of an optical cocoon tail behind the runaway HD 185806

Z. T. Spetsieri,1,2 P. Boumis,${}^{2\href
https://orcid.org/0000-0002-4260-940X{}}$ A. Chiotellis,${}^{2\href
https://orcid.org/0000-0002-4260-940X{}}$ S. Akras,${}^{2\href
https://orcid.org/0000-0003-1351-7204{}}$ S. Derlopa,2 S. Shetye,1 D. M.-A.
Meyer,3 D. M. Bowman,${}^{4\href https://orcid.org/0000-0001-7402-3852{}}$ V.
V. Gvaramadze5,6,7
1 Institute of Physics, Laboratory of Astrophysics, École Polytechnique
Fédérale de Lausanne (EPFL), Observatoire de Sauverny, 1290 Versoix,
Switzerland
2 Institute for Astronomy, Astrophysics, Space Applications and Remote
Sensing, National Observatory of Athens, 15236 Penteli, Greece
3 Universität Potsdam, Institut für Physik und Astronomie, Karl-Liebknecht-
Strasse 24/25, 14476 Potsdam, Germany
4 Institute of Astronomy, KU Leuven, Celestijnenlaan 200D, B-3001 Leuven,
Belgium
5 Sternberg Astronomical Institute, Lomonosov Moscow State University,
Universitetskij Pr. 13, Moscow 119992, Russia
6 Space Research Institute, Russian Academy of Sciences, Profsoyuznaya 84/32,
Moscow 117997, Russia
7 E. Kharadze Georgian National Astrophysical Observatory, Abastumani 0301,
Georgia
zoi.spetsieri<EMAIL_ADDRESS>to V.G. who passed away on 2
Sept. 2021

(Accepted 2022 June 29. Received 2022 June 24; in original form 2022 May 26)

###### Abstract

Studies on the circumstellar structures around evolved stars provide vital
information on the evolution of the parent star and the properties of the
local interstellar medium. In this work, we present the discovery and
characterization of an optical cocoon tail behind the star HD 185806.

The cocoon apex emission is puzzling, as it is detected in the infrared but
shows no signal in the optical wavelength. The H$\alpha$ and [O iii] fluxes of
the nebular structure vary from 2.7 to 8.5$\times$$10^{-12}$ erg s-1 cm-2 and
from 0.9 to 7.0$\times$$10^{-13}$ erg s-1 cm-2, respectively. Through high-
resolution spectroscopy, we derive the spectral type of the star, construct
the position-velocity diagrams of the cocoon tail for the H$\alpha$, [O iii]
and [N ii] emission lines, and determine its velocity in the range of $-$100
to 40 km s-1. Furthermore, we use SED fitting and MESA evolutionary models
adopting a distance of 900 pc, and classify HD 185806 as a 1.3 M⊙ star, in the
transition phase between the RGB and early AGB stages. Finally, we study the
morpho-kinematic structure of the cocoon tail using the astronomical software
SHAPE. An ellipsoidal structure, with an inclination of $\sim 19$° with
respect to the plane of sky is found to better reproduce the observed cocoon
tail of HD 185806.

###### keywords:

stars: RGB and early-AGB; stars: individual: HD 185806; stars: kinematics and
dynamics; ISM: kinematics and dynamics

††pubyear: 2022††pagerange: Discovery of an optical cocoon tail behind the
runaway HD 185806–LABEL:lastpage

##  1 Introduction

Several kinds of stars throughout the Hertzsprung-Russell (HR) diagram are
characterized by intense mass outflows in the form of stellar winds. These
winds drive gas and dust into the interstellar medium (ISM) and shape extended
circumstellar structures around the parent star (Langer, 2012; Smith, 2014).
Detailed studies of these circumstellar structures provide crucial insights on
the nature and evolution of the mass-losing star, the properties of the local
ISM and help us determine the impact of stellar outflows on the chemical and
dynamical evolution of their host galaxies.

Wind shaped circumstellar bubbles have been found to display a variety of
morphologies that deviate from the spherical symmetry being elliptical,
bipolar, multipolar and revealing arcs, lobes, blobs, jets, hydrodynamical
instabilities. Binary interactions (De Marco, 2009; van Marle et al., 2012;
Staff et al., 2016; Decin, 2019), magnetic field activity (Garcia-Segura,
1997; García-Segura et al., 2005), stellar rotation (Hoefner et al., 1996;
García-Segura et al., 1999; García-Segura et al., 2016) are only some of the
drivers of the morphologies formed around mass-losing stars. In addition,
facts related to the ambient medium such as density gradients/discontinuities,
magnetic fields, surrounding knots and clouds, pre-existing circumstellar
structures etc., further contribute on the final outcome of the resulted
circumstellar structures (e.g. Dwarkadas et al., 1996; Raga et al., 1998;
Wareing et al., 2007; Esquivel et al., 2010; van Marle et al., 2015).

In the case where the mass-losing star is moving supersonically, with respect
to the local ambient medium, a bow-shaped circumstellar structure is formed
resulting by the balance of the wind’s material and the ISM’s ram pressures
(Comeron & Kaper, 1998). Hence, bow-shaped wind structures are unique
laboratories for determining or at least constraining the mass loss history of
the central star, the surrounding ISM density and finally the stellar systemic
motion (e.g. Huthoff & Kaper, 2002; Gvaramadze et al., 2014; Meyer et al.,
2014; Chiotellis et al., 2016; Meyer et al., 2016, 2020).

Two main mechanisms have been proposed to be responsible for the formation of
circumstellar structures around supersonically moving stars whose space
velocities exceed the 30 $\rm km~{}s^{-1}$ a value that roughly corresponds to
three times the typical stellar space velocity (e.g. Blaauw, 1956; Gies &
Bolton, 1986; Huthoff & Kaper, 2002). The first mechanism implies close
dynamical interactions of either single stars (Poveda et al., 1967), binary-
binary interactions (Spitzer & Varshalovich, 1980) or finally by triple body
encounter between a binary and a single star (Gvaramadze & Gualandris, 2011)
that occur in a dense stellar environment such as in globular clusters. The
second mechanism suggests that the fast moving star is a member of a binary
system, which remained bound after the supernova explosion of the companion
star (Zwicky, 1957; Blaauw, 1961). These two mechanisms are more likely to
occur in massive stars. Indeed, the vast majority of bow shocks have been
observed in massive OB runaway stars. In these stars, the dust and gas are
accumulated behind the bow shock and are heated by the ultraviolet radiation
of the massive star (e.g. Henney & Arthur, 2019). Consequently, the dust
radiates at the infrared band (Draine & Lee, 1984), while the
shocked/photoionised gas produces optical line emission, such as H$\alpha$ and
[O iii]. For this reason, bow shocks are mostly detected in the mid-infrared
(e.g. van Buren et al., 1995; Noriega-Crespo et al., 1997; Peri et al., 2012;
Peri et al., 2015; Kobulnicky et al., 2016, 2017) and -less frequently- in the
optical band (e.g. Gull & Sofia, 1979; Kaper et al., 1997; Brown & Bomans,
2005; Boumis et al., 2009). In a few cases, non-thermal emission of bow shocks
has also been detected in the radio (syncrotron) and X-ray (inverse Compton)
bands (e.g. BD+43o3654, AE Aurigal; Benaglia et al., 2010; López-Santiago et
al., 2012, respectively).

In the low stellar mass regime, nebular structures are mostly observed around
stars in the red giant branch (RGB) and asymptotic giant branch (AGB) as
during this phase substantial mass outflows are met (e.g. Wenger et al., 2000;
Weigelt et al., 2002; Schöier et al., 2006; Castro-Carrizo et al., 2010;
Olofsson et al., 2010; Ali et al., 2013; Silva Aguirre et al., 2020). Most of
these bow shocks have been detected in the mid-infrared band. A large
Herschel/PACS survey of circumstellar structures around AGB stars revealed
that bow shocks are a common feature met in about 40$\%$ of the sample (Cox et
al., 2012). The stagnation point -i.e. nearest distance between the central
star and the circumstellar structure- was found to be smaller than 1 pc in all
cases. Such a result agrees with the mass loss history expected from AGB
stars, adopting typical Galactic warm neutral and ionised gas properties.

This paper focuses on the circumstellar structure around the M-type, long
period variable star (Kukarkin et al., 1970; Kato, 1999; Alfonso-Garzón et
al., 2012) at the constellation of Aquila, HD 185806 also known as V1279 Aql.
This star had been observed with the Wide-field Infrared Survey Explorer
(WISE; Wright et al., 2010) and displays a peculiar circumstellar structure
that resembles a "tail" feature. However no further investigation has been
conducted to identify the nature of the star or the properties of the
surrounding nebula.

In this work, we present the discovery of an optical cocoon tail behind of HD
185806. We obtain narrow band imaging in several filters and high resolution
optical spectroscopy, carrying out the first morpho-kinematic analysis of this
cocoon tail in the optical regime. The paper is structured as follows: In
Section 2, we present the optical images and high dispersion spectra of the
cocoon tail of HD 185806. The photometric and spectroscopic results of our
analysis are shown in Section 3. An extended discussion regarding the
properties of the star and its coccoon-tail, as well as, a 3D
morphokinematical model of the circumstellar structure are presented in
Section 4. We end with our conclusions in Section 5.

Table 1: Observing log. Columns from left to right specify the filter, the
central wavelength $\lambda_{\rm c}$, the wavelength difference
($\Delta\lambda$), the exposure time, the Observatory details and the date of
the observations. HIGH-RESOLUTION IMAGING with Aristarchos
---
Filter | $\lambda_{\rm c}$ | $\Delta\lambda$ | Exposure time | Observatory | Date |
| (Å) | (Å) | (s) |  |  |
H$\alpha$ \+ [N ii] 6584 Å | 6578 | 40 | 1800 | 2.3m Aristarchos | 2014 Aug 25-27 |
H$\alpha$ 6563 Å | 6567 | 17 | 1800 | 2.3m Aristarchos | 2014 Aug 25-27 |
[N ii] 6584 Å | 6588 | 17 | 1800 | 2.3m Aristarchos | 2014 Aug 25-27 |
[O iii] 5007 Å | 5011 | 30 | 1800 | 2.3m Aristarchos | 2014 Aug 25-27 |
HIGH-RESOLUTION SPECTROSCOPY with HERMES
V | 3800–9000 | - | 2000 | 1.2m Mercator Telescope | 2021 July 29 |
HIGH-RESOLUTION SPECTROSCOPY with MES
Area | Slit centres | Exp. time | Slit P.A./ width | Observatory | Date
| $\alpha$ | $\delta$ | (sec) | (degrees / $\mu$m) |  |
Slit 1 | 19h40m44.2s | 02°30′09″ | 1800 | 30 / 300 | 2.1m SPM | 2014 Nov 27
Slit 2 | 19h40m41.8s | 02°30′57″ | 1800 | 57 / 300 | 2.1m SPM | 2014 Nov 29

##  2 Observations

###  2.1 High-resolution imaging

Wide-field optical images of HD 185806 and its neighborhood were obtained with
the 2.3m Aristarchos telescope at Helmos Observatory in Greece in 2014. The
whole cocoon was imaged with four 5$\times$5 arcmin2 pointings. A CCD detector
with 1024$\times$1024 pixels2 with 24$\mu$m pixel size, resulting to 0.28
arcsec pixel-1 resolution, was used for the observations. Exposures of 1800 s
duration were obtained through the H$\alpha$, H$\alpha$+[N ii] and [N ii]
emission line filters on the main cocoon area, while the full mosaic was taken
only with the [O iii] filter (Boumis et al., 2016). In the latter, all fields
were projected on to a common origin on the sky and were subsequently combined
to create the final mosaic. The analysis of the images was carried out using
the standard IRAF routines for the bias subtraction, flat fielding correction,
bad pixel and column masking, sky background estimation and cosmic-rays
cleaning. The flux and astrometric calibration for all data frames was
performed using the STScI Guide Star Catalogue II (Lasker et al., 2008). We
performed photometry of all sources appearing in the images although we were
mainly focused on HD 185806. This was done in order to calibrate the magnitude
and flux of the star based on a standard zero point and have a larger sample
of stars to estimate the photometric errors. The photometric signal to noise
ratio (S/N) of our images was S/N$>$15 and the chi square value
$\chi{{}^{2}}\sim$2.0.

The resulting mosaic of the [O iii] image, where HD 185806 and the whole
cocoon tail can be seen is in Fig. 1, while in Fig. 2 we show the comparison
between the [O iii] and the H$\alpha$ emission line images. The H$\alpha$+[N
ii] and [N ii] images are not shown here, since the H$\alpha$+[N ii] and
H$\alpha$ images are almost identical, with only a very few small regions of
difference, where there is faint [N ii] emission.

###  2.2 High-dispersion and high resolution long-slit spectroscopy

High-dispersion long-slit echelle spectra were obtained at the 2.1 m telescope
in the Observatorio Astronómico of San Pedro Martir in Baja California,
Mexico, in its f/7.5 configuration using the Manchester Echelle Spectrometer
(Meaburn et al., 2003). The observations were carried out on 2014 November 27
and 29 (Table 1).

A 300 $\mu$m ($\equiv$ 3.9 arcsec and 20 km s-1) wide slit was used,
orientated at position angles of 30° and 57°. The 512 increments of the
2$\times$2 binned SITE CCD detector, each 0.624 arcsec long, give a total
projected slit length of 5.32 arcmin on the sky. In this spectroscopic mode
MES–SPM has no cross–dispersion consequently, for the present observations, a
filter of 90-Å bandwidth was used to isolate the 87th echelle order containing
the H$\alpha$ and [N ii] 6548, 6584 Å nebular emission lines, and a filter of
60-Å bandwidth to isolate the 114th echelle order containing the [O iii] line.
Integrations of 1800 s were obtained respectively for slit positions 1 and 2
marked in Fig. 3. The position–velocity (pv) arrays of H$\alpha$, [N ii] and
[O iii] line profiles from slit positions 1 and 2 are shown in Fig. 4.

We conducted similar pre-processing steps as for the images such as bias
subtraction, flat-fielding and illumination correction, cosmic-ray cleaning
and background subtraction. The background subtraction in our case was
conducted with extreme caution as over-subtraction of the background in
filters with fainter emission would lead to erroneous assumptions about the
spectral features of the bow shock structure. For the wavelength calibration
of the spectra the Th-Ar arc-lamp (Meaburn et al., 2003) was used.

We additionally collected a high resolution spectrum of HD 185806 with the
high-resolution fibre-fed HERMES spectrograph mounted on the 1.2-m Mercator
Telescope at the Roque de Los Muchachos Observatory, La Palma (Canary Islands;
Raskin et al. 2011). HERMES spectra have a spectral resolution of R = 85 000
and a wavelength coverage 380–900 nm. The spectrum was automatically reduced
through a dedicated pipeline at the end of the observing night. The
observations were obtained in the V band with a signal-to-noise ratio of 100
to ensure the accurate determination of the stellar parameters and of the
chemical abundances. A summary and log of our imaging and spectral
observations is provided in Table 1.

Figure 1: Mosaic image of the region surrounding HD 185806 in [O iii]
emission. Shadings run linearly from 0 to 8.5$\times$$10^{-12}$ erg s-1 cm-2.
Figure 2: A comparison between the [O iii] (left) and the H$\alpha$ (right)
images of the same region of HD 185806. The faint ring structures in the
H$\alpha$ image are ghosts due to the brightness of the star. Shadings run
linearly from 0 to 8.5$\times$$10^{-12}$ erg s-1 cm-2 and 0 to
7.5$\times$$10^{-13}$ erg s-1 cm-2, respectively. The images scale is
identical to that of Fig. 1. North is to the top, East to the left. Figure 3:
Full mosaic image of HD 185806 in [O iii] emission with slit positions 1 and 2
(in red). Shadings run linearly from 0 to 8.5$\times$$10^{-12}$ erg s-1 cm-2.
North is to the top, East to the left. Figure 4: Position-velocity images of
slit positions 1-2 in the emission light of H$\alpha$, [N ii] and [O iii].
Figure 5: Illustration of the quality of the match between observed and
synthetic spectra obtained for HD 185806 around the Zr line at 7849.37Å. The
upper panel presents a $\pm$ 3 Åzoom. The non-detection of Zr confirms that
the star does not show any s-process enrichment. Figure 6: The correlation
between the [O iii] and the W4 band at 22$\mu$m (yellow lines). The WISE
contours produced by a combined image from the three different wavebands,
where the star is bright in the first two and the nebula to the third. The
contours scale linearly from 1$\times 10^{-5}$ counts arcmin-2 s-1 to 3$\times
10^{-4}$ counts arcmin-2 s-1, with step 2.8$\times 10^{-5}$ counts arcmin-2
s-1.

##  3 Results

Table 2: Blue and red-shifted radial velocities in different positions across the slits, as shown in Fig. 3. In the last column are the minimum and maximum velocities measured in each slit for the different emission lines. Slit | Position | Emission line | Vhel | Vhel min-max
---|---|---|---|---
|  | ($\rm\AA$) | ( km sec-1) | ( km sec-1)
1 | A | H$\alpha$ |  $-75$ to $20$ |
| B |  |  $-100$ to $40$ |  $-100$ to $40$
| C |  |  $-85$ to $20$ |
1 | A | [N ii] 6584 |  $-50$ to $5$ |
| B |  |  |  $-75$ to $35$
| C |  |  $-65$ to $-7$ |
2 | D | H$\alpha$ |  $-55$ to $10$ |
| E |  |  $-85$ to $30$ |  $-100$ to $30$
| F |  |  $-90$ to $-20$ |
2 | D | [N ii] 6584 |  $-50$ to $5$ |
| E |  |  |  $-80$ to $5$
| F |  |  $-80$ to $-30$ |
2 | D | [O iii] 5007 |  $-55$ to $15$ |
| E |  |  $-80$ to $-40$ |  $-90$ to $15$
| F |  |  $-60$ to $-20$ |

###  3.1 The optical emission line images

Figure 1 displays the optical [O iii] image of HD 185806 surrounded by an
elliptical shape nebula, with a filamentary structure developed. The brighter
filaments have [O iii] fluxes from 7.1 to 8.5$\times$$10^{-12}$ erg s-1 cm-2
and the fainter ones from 2.7 to 3.1$\times$$10^{-12}$ erg s-1 cm-2. For
comparison, the [O iii] and H$\alpha$ emission of the nebula is shown in Fig.
2. The nebula is evident in H$\alpha$ displaying the same structure as in [O
iii], however all filaments are fainter at this wavelength. The bright
filaments have H$\alpha$ fluxes between 3.5 and 7.0$\times$$10^{-13}$ erg s-1
cm-2 and the fainter ones between 8.5$\times$$10^{-14}$ erg s-1 cm-2 and
1.5$\times$$10^{-13}$ erg s-1 cm-2. The [O iii]/H$\alpha$ flux ratio of the
nebula is in the range from 10 to 30 being higher in the fainter filaments.

It is interesting to note that the emission is stronger in the opposite outer
sides of the cocoon (areas A, C and D of Fig. 3 ) between $\alpha\simeq$
19h40m45s, $\delta\simeq$ 02°31′30″; $\alpha\simeq$ 19h40m47s, $\delta\simeq$
02°31′10″ and $\alpha\simeq$ 19h40m40s, $\delta\simeq$ 02°30′00″;
$\alpha\simeq$ 19h40m44s, $\delta\simeq$ 02°29′25″. There are also two bright
filaments almost perpendicular to the outer ones at $\alpha\simeq$ 19h40m46s,
$\delta\simeq$ 02°30′15″ (close to area B) and around $\alpha\simeq$
19h40m42s, $\delta\simeq$ 02°30′30″ (close to area E). On the contrary, both
the bow-shock and the cocoon tail regions opposite to the star have either
faint or no obvious emission.

###  3.2 High-resolution spectra of the optical cocoon tail around HD 185806

The full mosaic of HD 185806 with the cocoon tail surrounding it, is shown in
Fig. 3. The slit positions are indicated in red color while the regions where
the heliocentric velocity of the gas was derived (Table 2) are labeled as A,
B, C, D, E, and F. Figure 4 shows the position velocity (PV) diagrams in the
corresponding emission lines for slit 1 (S1) (top two panels) and slit 2 (S2)
(lower three panels). A general result is that in all positions and in both
slits, the whole structure of the cocoon reveals a bulk velocity towards us.
Furthermore, at the same emission lines, there is a good agreement in
velocities in positions from the same or nearby filaments, as it is for A and
D, B and E, C and F. In Table 2, we also list the minimum and maximum
velocities measured in each emission line, where it can be seen that the
H$\alpha$ ranges from $-$100 to 40 km s-1, the [N ii] from $-$80 to 35 km s-1
and the [O iii] from $-$90 to 15 km s-1.

It is also noticeable in the PV diagrams that the H$\alpha$ and [O iii]
emissions are detected all over the cocoon area, while in the [N ii] emission
are detected only some bright spots (A/C and D/F for the S1 and S2 slit
positions). The detection of the H$\alpha$, [N ii], and [O iii] emission lines
in the cocoon tail of HD 185806 in conjunction with the low temperature of the
star implies the presence of a shock-heated gas. If the orientation of the
cocoon is taken into account, the expansion velocity of the ionized gas would
be $>$ 70-80 km s-1 which can explain the strong [O iii] emission. However,
follow-up low-dispersion spectroscopic data are required in order to provide
insights on the physical conditions of the gas.

###  3.3 High-resolution spectrum of HD 185806

The high resolution spectrum from HERMES yielded that HD 185806 is an M4-type,
long period variable star. We estimated the atomic abundances following the
same method as described in Shetye et al. (2018). The spectra of S and M stars
are dominated by molecular bands making the detection of unblended atomic
s-process lines quite rare (Smith & Lambert, 1990). This is why a spectral-
synthesis approach is required, as opposed to relying solely on equivalent
widths. In order to check whether the star is undergoing active slow neutron
capture nucleosynthesis (s-process), which is typical in thermally-pulsing AGB
stars, we checked for any presence of Technetium (Tc). Tc is an s-process
element with no stable isotope. One of the isotopes 99Tc of the star has a
half-life of 210 000 yrs. This peculiar property of Tc makes it the best
indicator for any active s-process nucleosynthesis in S and M type stars
(Shetye et al., 2018). We checked the three Tc I lines located at 4239.191 Å,
4267.27 Å and 4297.06 Å, and did not detect any absorption of Tc. We conclude
that HD 185806 is not a thermally-pulsing AGB star. As in Shetye et al.
(2018), we used the two Zr lines at 7819.37 Å and 7849.37 Å with transition
probabilities from laboratory measurements (Biemont et al., 1981). In the case
of HD 185806 we do not report any overabundance of Zr, confirming that the
star does not show any enrichment of s-process elements. We derived the [Fe/H]
using some Fe lines at [Fe/H]= -0.50 dex revealing the sub-solar metallicity
of the star. We additionally derived the nitrogen abundance from the CN lines
in the 8000-8100 Årange. In particular, using the lines listed in Merle et al.
(2016). We report a nitrogen abundance [N/Fe] = +0.6 to +1 dex. An
illustration of the observed and synthetic spectra obtained for HD 185806
around the Zr line at 7849.37 Å are shown in Fig. 5.

Table 3: Stellar properties of HD 185806. RA | Dec | Teff | Luminosity | mass | parallax | radial velocity | transverse velocity | systemic velocity
---|---|---|---|---|---|---|---|---
(J2000) | (J2000) | K | L⊙ | M⊙ | mas | km s-1 | km s-1 | km s-1
19:40:47.56 | +02:29:28.64 | 3400$\pm$100 | 1548 $\pm$200 | 1.3 | 1.109 | 14.0 | 43.75 | 45.0

###  3.4 WISE observations

HD 185806 and its cocoon tail have also been observed with WISE. The stellar
source is detected in all bands (W1 (3.4 $\mu$m), W2 (4.6 $\mu$m), W3 (12.1
$\mu$m), and W4 (22 $\mu$m) (Wright et al., 2010) but it is saturated in the
first two. The cocoon tail is detected only in the W4 band which indicates the
presence of warm dust. Fig. 6 shows the [O iii] image with the contours from
the W4 image over plotted in yellow. The contours surround both the star and
the nebula, as well as, the interior filaments of the nebula revealing a
strong correlation between the optical and infrared emission. This correlation
is not uncommon, since it has been found in other cases, such as in Wolf-Rayet
star nebulae (Toalá et al., 2015). The W4 image shows that the emission in
this band is dominated by thermal continuum emission from dust spatially
coincident with the thin nebular shell or most likely with the leading edge of
the nebula. Intriguingly, the bow shock that is absent from the optical images
of HD 185806, is likely detected in the W4 band as it can be seen from the
outer contours (Fig. 6).

This optical/infrared discrepancy at the location of the forward shock of the
cocoon may partly originate from the penetration of interstellar dust into the
inner circumstellar medium of HD 185806. This engenders a shift between the
location of the outer layer of the cocoon and the location where the dust
accumulates, as shown in van Marle et al. (2011) in the context of the cool
runaway red supergiant star Betelgeuse. Finally, another deviation between the
optical and infrared occurs at the edge of the cocoon tail, opposite to the
star, where the faint [O iii] filament is not accompanied by IR emission.

##  4 Discussion

###  4.1 Properties and nature of HD 185806

Our spectroscopic analysis confirmed that HD 185806 is an M4-type long period
variable star. Based on Gaia eDR3 we estimated the distance of HD 185806 at
D=909 pc (Gaia Collaboration et al., 2016, 2018). The angular velocity of the
star is observed to be $\mu\sim 10.2~{}\rm mas~{}yr^{-1}$ and for the given
distance corresponds to a transverse stellar velocity of $V_{T}=43.75~{}\rm
km~{}s^{-1}$. Finally, the Gaia eDR3 radial velocity value is found to be
$V_{r}=14.0~{}\rm km~{}s^{-1}$.

We used the photometric data available in virtual observatory spectral energy
distribution (SED) analyzer (VOSA; Bayo et al., 2008) to determine the
effective temperature (Teff) and luminosity (L) of the star. VOSA offers tools
for the SED analysis, allowing the estimation of the stellar parameters. The
synthetic photometry used by VOSA is calculated by convolving the response
curve of the used filter set with the theoretical synthetic spectra. Then, a
statistical test is performed, via $\chi{{}^{2}}$ minimization, to estimate
which set of synthetic photometry best reproduces the observed data. The
stellar parameters are obtained from the comparison of the observed spectrum
with a simulated one with known stellar parameters. In our case, using VOSA we
created the SED for HD 185806 and fitted the data with a library of high-
resolution synthetic stellar spectra (Coelho et al., 2005) created by the
PFANT code (Cayrel et al., 1991; Barbuy et al., 2018). The SED generated for a
distance of 900 pc is shown in Fig. 7 and corresponds to Teff= 3400$\pm$100 K
and L=1548$\pm$200 L⊙. The parameters derived from VOSA seem to be in well
agreement with the spectra. The extracted properties of HD 185806 are
summarized in Table 3.

We employed the Modules for Experiments in Stellar Astrophysics (MESA)
Isochrones and Stellar Tracks (MIST; Choi et al., 2016) to estimate the mass,
age, and evolutionary state of HD 185806 considering solar metallicity. MESA
models use solar-scaled abundances covering a wide range of ages, masses,
phases, and metallicities computed within a single framework. The models were
generated assuming rotation at 40 per cent of the critical breakup speed.
Based on the spectral type, effective temperature, luminosity and [Fe/H], we
plotted evolutionary tracks of 1-3 M⊙ with steps of 0.1 M⊙, using a Galactic
extinction of AV=1.12 mag by NASA/IPAC extra-galactic database (NED; Schlegel
et al., 1998). The evolutionary track which best fits the stellar parameters
of HD 185806 corresponds to an 1.3 M⊙ star. Figure 8 illustrates the
evolutionary track of a 1.3 M⊙ star, in the HR diagram, together with the
extracted parameters of HD 185806. The star is lying on the red giant branch
(RGB) phase. This is consistent with the chemical signatures of the star that
were discussed in Section 3. Such an evidence is also aligned to the formation
of the cocoon tail detected around HD 185806 as RGB stars are characterised by
intense and slow stellar winds, able to shape dense circumstellar structures
around the star.

Figure 7: Spectral Energy Distribution derived by VOSA. Colour-coded are the
various photometric points from different surveys. Observations from 2MASS
(Skrutskie et al., 2006) are shown in green, from IRAS (Neugebauer et al.,
1984) in orange, from the AKARI (Ishihara et al., 2010) in yellow, from WISE
(Wright et al., 2010) in pink, from Tycho (Høg et al., 2000) in grey, from
APASS (Henden et al., 2016) in purple, from Gaia eDR3 (Gaia Collaboration et
al., 2016, 2018) in magenta, from Gaia eDR3 (Gaia Collaboration et al., 2021)
in brown andhttps://www.overleaf.com/project/6048838257fb47049f965521 Pan-
Starrs PS1 (Chambers & et al., 2017) in light green. The red curve is the SED
fitting generated by high-resolution synthetic stellar spectra by (Coelho et
al., 2005)

.

Figure 8: Hertzsprung-Russell diagram generated from MESA for a 1.3 M⊙ star.
The phases of stellar evolution are displayed in different colors. The main
sequence is shown in blue color, the red giant branch (RGB) is in orange, the
Core He burning phase is shown in green, the early AGB phase in yellow, the
thermal pulsating AGB phase is in red and the post AGB phase is in black
color. HD 185806 is shown as a red star. Figure 9: Mass-loss vs. r0 plot. The
four panels show the distance of the bow shock from the star for varying
values of wind velocity (Uwind) and ISM density (nism). The blue line
corresponds to a nism=0.1, the orange line to nism=1 and the green line to
nism=10. The dashed red line shows the distance below which the bow shock is
not observable.

###  4.2 On the morphology of the cocoon tail around HD 185806

The evolutionary state and the high spatial velocity of HD 185806 implies that
the observed wind-blown bubble morphology around the star is dominantly
determined by the interaction between the stellar wind material and the
surrounding ISM.

The circumstellar structures around supersonically moving stars are
characterized by a bow-shaped shell formed in the direction of the stellar
motion, followed by a cocoon/tail morphology. In the case where the wind and
the ISM flows are supersonic, with respect to the sound speed of the local
medium, shock waves are formed. The formed shock waves compress and heat up
the medium that are propagating. The resulting bow shock structure consists of
a forward shock separating the unshocked and shocked ISM, a wind termination
shock that marks the border between the free-expanding wind and the shocked
wind, and finally a contact discontinuity lies between the two shocks
separating the shocked wind from the shocked ISM shells (see Weaver et al.
1977 and Lamers & Cassinelli 1999). The position and shape of the bow shock is
governed by the ram pressure balance between the stellar wind and the ISM,
that in turn is determined by the wind mass loss rate (${\bm{\dot{}}{M}}$) and
terminal velocity ($u_{w}$), the stellar systemic velocity ($u_{*}$) and the
local ISM density ($\rho_{\rm ISM}$). In the direction of the stellar motion
the bow shock is at its shorter distance - called as stagnation point - the
radius of which is given by (Baranov et al., 1971; Dyson, 1975; Wilkin &
Stahler, 1996):

$R_{0}=\sqrt{\frac{{\bm{\dot{}}{M}}\ v_{w}}{4\pi\ \rho_{ism}\ v^{2}_{\star}}}$
(1)

The CSM around HD 185806 displays a cocoon-tail structure behind the star
which indicates the presence of shocked ISM/wind material by the forward and
termination shock of the wind bubble. However, even though a bow shock
structure is clearly depicted at the WISE images (Wright et al., 2010) (see
the yellow contours in Fig. 6), there is no evidence of the bow shock in the
optical band. The current complex properties of the CSM around HD 185806 and
the absence of an optically bright bow shock shell can be attributed to the
following three main reasons (or to a combination of them):

####  4.2.1 Proximity of the bow shock to the central star

In cases where the formed bow shock lies very close to the star it is likely
to be out-shined by the stellar radiation or/and to be obscured by the
produced stellar dust. To assess the plausibility of this statement in Fig. 9
we illustrate the radial distance between the star and the stagnation point,
as given by Eq. (1), for different values of RGB stellar winds and of the ISM
density. The dashed red line in Fig. 9 indicates the threshold below which the
bow shock is not expected to be observed in the optical as it coincides with
the area where the central star radiation dominates. As it is clearly
portrayed, for the case where the ISM density is higher $n_{\rm ISM}\geq
1cm^{-3}$ the predicted stagnation point radius is below the defined threshold
for moderate mass loss rates ($\dot{M}\lesssim 10^{-7}~{}\rm
M_{\odot}~{}yr^{-1}$), something that makes this scenario plausible.

From a statistical point of view, half of the RGB/AGB stars and red
supergiants with bow-shock shells studied by Cox et al. (2012) have a
stagnation point at a de-projected distance of $<$0.08 pc, while there are 29
more stars without any detection. Therefore, the non-detection of a clear bow-
shock in the HD 185806 is not unique or even rare, something that advocates
further towards the credibility of the initial assumption.

####  4.2.2 Time variable stellar wind properties

A second possible explanation of the non-detection of an optical emitting bow
shock is based on the current evolutionary state of HD 185806. As we discussed
in Section 4.1, the observed properties of the central star indicate that HD
185806 is currently at the RGB phase evolving towards the core helium
burning/early AGB phase. If indeed the central star is in a transition state
between two evolutionary stages it is likely that in the recent past its
stellar wind properties (mass loss rate and wind velocity) have changed
substantially. In particular, in this transition phase the mass loss rate is
expected to decrease and the wind terminal velocity to increase. As a result,
the freely expanding wind density drops significantly. In the framework of
this scenario, the lack of an optically bright bow shock around HD 185806 is
attributed to the low wind-bubble density in which the terminating shock is
propagating. In addition, the information about the stellar wind properties
alteration is not communicated spontaneously throughout the circumstellar
structure but it travels within the dynamic timescale of the wind flow:
$t_{\rm flow}\sim\frac{r}{u_{w}}$, where $r$ the distance from the mass-losing
star and $u_{w}$ the wind velocity. Thus, there is the chance that we are
currently witnessing a transition phase of the CSM around HD 185806, were in
the region close to the stagnation point a low density wind is dominant, while
the observed cocoon tail structure corresponds the imprint of an anterior
phase where dense and bright CSM structures were formed.

####  4.2.3 Deviation of the ambient medium from homogeneity

The above discussed arguments assume an uniform ambient medium in the warm
phase ($T_{\rm ISM}\approx 1000-8000\,\rm K$) of the ISM of the Milky Way
(McKee, 1995). The establishment of a stellar wind bow shock being both a
function of the wind and ISM density combined with the absence of an observed
arc around runaway stars, can indicate runaway stellar objects moving through
a hot low-density region (huthoff_aa_383_2002). However, this argument does
not apply to HD 185806 since the bow shock is visible in near-infrared but not
in optical.

Nevertheless, any deviation from uniformity in the ambient medium affects the
development of its circumstellar nebula. For example, the presence of native
ISM overdensities and/or disturbances of the velocity field interacting with
the forward shock of the bow shock induces distortions of the layers of
shocked material that produce bright optical, free-free and near-infrared
knots in the regions of shocked material (Baalmann et al., 2021; Herbst et
al., 2022). Consequently, potential irregularities in the ISM density field in
the vicinity of HD 185806 can only take the form of local underdensities which
would make its bow shock deviate from the analytic solution of Wilkin (1996).
Such an example, is the production of a local extension of the forward shock,
which is affected by Rayleigh-Taylor and other non-linear instabilities
(Vishniac, 1994). The star could be moving in an ISM gradually decreasing
along the direction of stellar motion. This may explain the wise $22\,\mu\rm
m$ protuberances in the south-west region of Fig. 6.

Another assumption that has been tacitly made is that the ambient medium
around HD 185806 is characterized by the local ISM properties. Thus, it has
been neglected any ambient medium modifications that occurred in a prior phase
of the central star’s evolution, namely the main-sequence (MS) phase.

The complex internal morphology of the cocoon around HD 185806 could
potentially be explained considering that the current circumstellar structure
around HD 185806 resulted by the evolution of the RSG wind bubble within the
interior of the previously formed MS bubble. The arcs and filamentary shapes
visible via the [O iii] emission of Fig. 6 not only trace the internal border
of the cocoon, but also appear as bright structures parallel to the direction
of stellar motion. The former can reasonably be associated to the contact
discontinuity of the main-sequence wind-blown bubble at the moment it is
shocked by the new-born RGB wind. Regarding the twisted optical arcs
encompassed by the infrared cocoon we propose that they are projection effects
similar to those found in the three-dimensional magneto-hydrodynamical
simulations of Meyer et al. (2021). It is shown that the instabilities
developing in the termination shock of the stellar wind bow shock of red
supergiant stars can turn, by projection effects, into a large variety of
optical H$\alpha$ filaments, including parallel optical segments to the
direction of stellar motion, see Fig. 10f,h and Fig. 11e,h in Meyer et al.
(2021). We argue that a similar mechanism is at work inside of the cocoon of
HD 185806.

HD 185806 should have recently entered into its current, red giant branch,
evolutionary phase, its circumstellar medium is not in a steady state (Mohamed
et al., 2012) and the time-variability of its stellar wind is a key parameter
into the deep understanding of its cocoon. With this in mind it raises the
question of its future evolution, particularly within the predictive models
for fast-moving stars in the AGB, thermally-pulsating AGB stars and planetary
nebula phase of Villaver & Stanghellini (2005). The series of papers of
Villaver & Stanghellini (2005) strengthens the interpretation of the observed
cocoon as a transient circumstellar structure generated by the onset of the
RGB evolutionary phase, which might be replaced by an AGB wind to be released
into the cocoon. The surroundings of HD 185806 may then eventually evolve to
those of a thermal-pulsating AGB star, as depicted in the numerical models of
Villaver et al. (2012).

Figure 10: (a) The 3D morpho-kinematic (MK) model of HD overlaid on its
observational wide field image in [O iii]. The two black lines represent the
slits positions, while the red circle is the location of the star. (b) The 3D
MK model in grid representation. The red dot is the location of the star,
while the vectors are the velocity vectors that deduced from the velocity
field of the shell of the cocoon tail. For description, see in the subsection
4.4. (c) The observational PV diagrams in black deduced from the echelle
spectra in H$\alpha$ and [O iii] emission lines, overlaid with the synthetic
PV diagrams reproduced with SHAPE. Blue and red colours corresponds to the
blue- and red-shifted gas.

###  4.3 Properties of the cocoon’s forward shock

The exact understanding of the morphology of the infrared/optical forward
shocks of the cocoon around HD 185806 reveals a dynamical and internal
synthesis of a nebula which exceeds that of a simple steady-state stellar wind
bow shock. Further investigation would require more sophisticated numerical
models. However, we can still elaborate on the shock properties by means of
analytic arguments, as done in, e.g. Stevens et al. (1992); Gvaramadze et al.
(2011); Chiotellis et al. (2012, 2016); Meyer et al. (2017). In order to
characterise the observed shock, we will use the methodology described in
Gvaramadze et al. (2011), and the latter works which were designed to
characterise simulated hydrodynamical shocks.

Bow shocks are mainly characterised by two quantities, namely the dynamical
$t_{\rm dyn}=r_{0}/v$ and cooling $t_{\rm cool}=3k_{\rm B}T_{\rm ps}/n_{\rm
ps}\Lambda(T_{\rm ps})$ timescales, respectively, where $k_{\rm B}$ is the
Boltzmann constant, $T_{\rm ps}$ the post-shock temperature at the forward
shock, $n_{\rm ps}$ the post-shock density at the forward shock and
$\Lambda(T_{\rm ps})$ the corresponding cooling rate by optically-thin
processes, see Meyer et al. (2014). According to Fig. 9, the stand-off
distance of the cocoon is $\approx 0.1\leq r_{0}\leq 1.0\,\rm pc$. Since the
bulk motion of the star is $v_{\star}=45\,\rm km\,\rm s^{-1}$, we have
$v=v_{\star}/4\approx 11.25\,\rm km\,\rm s^{-1}$ and therefore $t_{\rm
dyn}\approx 0.0088$$-$$0.088\,\rm Myr$ assuming $r_{0}=0.1$$-$$1.0\,\rm pc$.

Similarly, the cooling timescale is obtained by estimating first the
conditions in the post-shock region at the forward shock using the relations
in section 5.3 of Gvaramadze et al. (2011). This gives $T_{\rm ps}\approx
1.4\,\times 10^{5}(v_{\star}/100\,\rm km\,\rm s^{-1})^{2}\approx 28350\,\rm K$
and $n_{\rm ps}=4n_{\rm ISM}=0.4$$-$$40\,\rm cm^{-3}$ (assuming the explored
parameter space of Fig. 9). Using the cooling curve detailed in Meyer et al.
(2014), one can find $\Lambda(T_{\rm ps})=4.78\times 10^{-23}\,\rm erg\,\rm
cm^{3}\,\rm s^{-1}$. Consequently, $t_{\rm cool}=0.019$$-$$0.00019\,\rm Myr$
for $n_{\rm ISM}=0.1$$-$$10\,\rm cm^{-3}$.

The creation of large bow shock ($r_{0}=1.0\,\rm pc$) is possible with a
diluted ambient only ($n_{\rm ISM}=0.1\,\rm cm^{-3}$), for which it is found
that $t_{\rm cool}=0.019\sim t_{\rm dyn}\approx 0.088\,\rm Myr$, indicating
that the shock is adiabatic, with both timescales being of the same order of
magnitude. The creation of small bow shock ($r_{0}=0.1\,\rm pc$) is possible
in a dense ISM ($n_{\rm ISM}=10\,\rm cm^{-3}$), gives $t_{\rm cool}=0.00019\ll
t_{\rm dyn}\approx 0.0088\,\rm Myr$, meaning that the shock is radiative.

The stellar evolution calculations performed in the context of HD 185806 show
that its mass-loss rate is $\dot{M}\leq 10^{-7}\,\rm M_{\odot}\,\rm yr^{-1}$
This value favors the scenario of the existence of a small bow shock
($r_{0}=0.1\,\rm pc$) and therefore of a radiative forward shock of the
circumstellar cocoon of HD 185806. This should additionally be confirmed by
multi-dimensional hydrodynamical calculations in our following studies.

###  4.4 3D morpho-kinematic modeling of the cocoon shell

For the Morpho-kinematic (MK) modeling of the cocoon tail of HD 185806 in
three-dimensions, we employed the astronomical software SHAPE (Steffen et al.,
2011). This code has been widely used to model morpho-kinematic structures of
planetary nebulae (e.g. Akras & Steffen, 2012; Akras & López, 2012; Akras et
al., 2015; Clyne et al., 2015), nova shells (Harvey et al., 2020) and more
recently complex supernova remnants (Derlopa et al., 2020). Despite the cocoon
tail of HD 185806 displays a simple ellipsoidal morphology, the intriguing
formation of the nebula behind the highly moving star motivated us to
reproduce the observed position-velocity diagram and reconstruct the 3D
velocity field.

The wide field image of HD 185806 in the [O iii] line, along with the PV
diagrams in the H$\alpha$ and [O iii] emission lines were used to constrain
the morpho-kinematic structure of the cocoon tail. An ellipsoidal structure
with a projected size of 3.6$\times$2.4 arcmin2 can adequately reproduce the
nebula behind the HD 185806 star. The position angle of the cocoon tail is
found to be $\sim$19° and the inclination angle $\sim$10°, both with respect
to the plane of the sky, i.e. the star goes inwards the page, as it is shown
in Fig. 10(a).

For the expansion velocity law of the shell gas, a homologous expansion law
($\overrightarrow{V}$ $=$ $k$ $\overrightarrow{r}$) was applied. The centre of
the velocity law is on the HD 185806 star and not at the geometrical centre of
the nebula such as in PNe, resulting in faster expanding gas at larger
distances from the star. Fig. 10(b) shows the 3D grid structure of the cocoon
tail along with the vectors of the velocity field (magnitude and direction).
The red dot represents the location of the star. Fig. 10(c) illustrates the
matching between the observed (black) and modelled (colored) PV diagrams
reproduced with SHAPE. According to our model, the mean value for the $k$
constant parameter in the homologous expansion is $31.3$ km sec-1 arcmin-1.

##  5 Conclusions

We have investigated the properties of the source HD 185806 and of its newly
detected optical cocoon tail. By performing a photometric and spectroscopic
analysis of the source in the optical filters H$\alpha$, [O iii], H$\alpha$+[N
ii] and [N ii] we report that:

  * •

The star is an M type star followed by an optical cocoon tail with intense [O
iii] emission. The [O iii] flux of the filaments ranges from 2.7 to
8.5$\times$$10^{-12}$ erg s-1 cm-2. The spectroscopic analysis of the results
shows that the heliocentric expansion velocity of the filaments ranges between
-100 to 40 km s-1, and the systemic velocity is
45 km s-1.

  * •

We used the virtual observatory SED analyzer (VOSA; Bayo et al., 2008) in
which we imported all the available photometric data of the star and
constructed the SED of the source from 4000 Å up to 100 $\micron$. The SED
implies an effective temperature Teff=3400$\pm$100 K and logL=3.18 L⊙ as well
as an infrared excess. The stellar parameters obtained from VOSA were used to
create an HR diagram using the MESA models and evolutionary tracks. MESA
models implied that the star is a 1.3 M⊙ star in the RGB/early-AGB phase.

  * •

We investigated the absence of an optically bright bow shock shell around HD
185806 which can be attributed to two main factors: the distance of the
stagnation point regarding the star and the current evolutionary state of HD
185806. In the first case the bow shock could be outshined by the stellar
radiation or/and be obscured by the large amount of dust that lies close to
the central star (Fig. 9). In the second case, as the models imply that the
star is in the transition phase between the RGB/early AGB phase, the wind
mass-loss rate undergoes a decrease and the wind terminal velocity increases.
As a result, the freely expanding wind density props by several orders of
magnitude and the dense stellar wind bubble is transformed into a low density
cavity. Regarding the statistical facts provided in Cox et al. (2012) the non-
detectability of a bow shock in the optical regime is not uncommon.

  * •

The 3D morpho-kinematic model we created for the cocoon tail shows that its
morphology can be reproduced by an ellipsoidal structure with a projected size
of 3.6$\times$2.4 arcmin2, a position angle of $\sim$19° and an inclination of
$\sim$10° (Fig. 10a). We studied the expansion velocity using a homologous law
($\overrightarrow{V}$ $=$ $k$ $\overrightarrow{r}$) as shown in Fig. 10(b)
and, according to our model, the mean value for the $k$ constant parameter is
$31.3$ km sec-1 arcmin-1.

## Acknowledgements

This research is co-financed by Greece and the European Union (European Social
Fund-ESF) through the Operational Programme “Human Resources Development,
Education and Lifelong Learning 2014-2020” in the context of the project “On
the interaction of Type Ia Supernovae with Planetary Nebulae” (MIS 5049922).
A.C. acknowledge the support of this work by the project “PROTEAS II” (MIS
5002515),which is implemented under the Action “Reinforcement of the Research
and Innovation Infrastructure”, funded by the Operational Programme
“Competitiveness, Entrepreneur- ship and Innovation” (NSRF 2014–2020) and co-
financed by Greece and the European Union (European Regional Development
Fund). SA acknowledges support under the grant 5077 financed by IAASARS/NOA.
DMB gratefully acknowledges a senior postdoctoral fellowship from the Research
Foundation Flanders (FWO) with grant agreement no. 1286521N. This work is
based on observations made with the “Aristarchos” telescope operated on the
Helmos Observatory by the Institute of Astronomy, Astrophysics, Space
Applications and Remote Sensing of the National Observatory of Athens. This
work has made use of data from the European Space Agency (ESA) mission Gaia
(https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and
Analysis Consortium (DPAC,
https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has
been provided by national institutions, in particular the institutions
participating in the Gaia Multilateral Agreement. Data products from the Wide-
field Infrared Survey Explorer, which is a joint project of the University of
California, Los Angeles, and the Jet Propulsion Laboratory/California
Institute of Technology, funded by the National Aeronautics and Space
Administration. The authors acknowledge the North-German Supercomputing
Alliance (HLRN) for providing HPC resources that have contributed to the
research results reported in this paper. The authors would like to dedicate
this article to Vasilii V. Gvaramadze who passed away on 2 September 2021,
since we started this project together.

## Data Availability

The data underlying this article will be shared on reasonable request to the
corresponding author.

## References

  * Akras & López (2012) Akras S., López J. A., 2012, MNRAS, 425, 2197
  * Akras & Steffen (2012) Akras S., Steffen W., 2012, MNRAS, 423, 925
  * Akras et al. (2015) Akras S., Boumis P., Meaburn J., Alikakos J., López J. A., Gonçalves D. R., 2015, MNRAS, 452, 2911
  * Alfonso-Garzón et al. (2012) Alfonso-Garzón J., Domingo A., Mas-Hesse J. M., Giménez A., 2012, A&A, 548, A79
  * Ali et al. (2013) Ali A., Ismail H. A., Snaid S., Sabin L., 2013, A&A, 558, A93
  * Baalmann et al. (2021) Baalmann L. R., Scherer K., Kleimann J., Fichtner H., Bomans D. J., Weis K., 2021, A&A, 650, A36
  * Baranov et al. (1971) Baranov V. B., Krasnobaev K. V., Kulikovskii A. G., 1971, Soviet Physics Doklady, 15, 791
  * Barbuy et al. (2018) Barbuy B., Trevisan J., de Almeida A., 2018, PFANT: Stellar spectral synthesis code (ascl:1812.003)
  * Bayo et al. (2008) Bayo A., Rodrigo C., Barrado Y Navascués D., Solano E., Gutiérrez R., Morales-Calderón M., Allard F., 2008, A&A, 492, 277
  * Benaglia et al. (2010) Benaglia P., Romero G. E., Martí J., Peri C. S., Araudo A. T., 2010, A&A, 517, L10
  * Biemont et al. (1981) Biemont E., Grevesse N., Hannaford P., Lowe R. M., 1981, ApJ, 248, 867
  * Blaauw (1956) Blaauw A., 1956, PASP, 68, 495
  * Blaauw (1961) Blaauw A., 1961, Bull. Astron. Inst. Netherlands, 15, 265
  * Boumis et al. (2009) Boumis P., Meaburn J., Lloyd M., Akras S., 2009, MNRAS, 396, 1186
  * Boumis et al. (2016) Boumis P., Akras S., Leonidaki I., Chiotellis A., Kopsacheili M., Alikakos J., Nanouris N., Mavromatakis F., 2016, in Supernova Remnants: An Odyssey in Space after Stellar Death. p. 15
  * Brown & Bomans (2005) Brown D., Bomans D. J., 2005, A&A, 439, 183
  * Castro-Carrizo et al. (2010) Castro-Carrizo A., et al., 2010, A&A, 523, A59
  * Cayrel et al. (1991) Cayrel R., Perrin M. N., Buser R., Barbuy B., Coupry M. F., 1991, A&A, 247, 122
  * Chambers & et al. (2017) Chambers K. C., et al. 2017, VizieR Online Data Catalog, p. II/349
  * Chiotellis et al. (2012) Chiotellis A., Schure K. M., Vink J., 2012, A&A, 537, A139
  * Chiotellis et al. (2016) Chiotellis A., Boumis P., Nanouris N., Meaburn J., Dimitriadis G., 2016, MNRAS, 457, 9
  * Choi et al. (2016) Choi J., Dotter A., Conroy C., Cantiello M., Paxton B., Johnson B. D., 2016, ApJ, 823, 102
  * Clyne et al. (2015) Clyne N., Akras S., Steffen W., Redman M. P., Gonçalves D. R., Harvey E., 2015, A&A, 582, A60
  * Coelho et al. (2005) Coelho P., Barbuy B., Meléndez J., Schiavon R. P., Castilho B. V., 2005, A&A, 443, 735
  * Comeron & Kaper (1998) Comeron F., Kaper L., 1998, A&A, 338, 273
  * Cox et al. (2012) Cox N. L. J., et al., 2012, A&A, 537, A35
  * De Marco (2009) De Marco O., 2009, PASP, 121, 316
  * Decin (2019) Decin L., 2019, in ALMA2019: Science Results and Cross-Facility Synergies. p. 53, doi:10.5281/zenodo.3585272
  * Derlopa et al. (2020) Derlopa S., Boumis P., Chiotellis A., Steffen W., Akras S., 2020, MNRAS, 499, 5410
  * Draine & Lee (1984) Draine B. T., Lee H. M., 1984, ApJ, 285, 89
  * Dwarkadas et al. (1996) Dwarkadas V. V., Chevalier R. A., Blondin J. M., 1996, ApJ, 457, 773
  * Dyson (1975) Dyson J. E., 1975, Ap&SS, 35, 299
  * Esquivel et al. (2010) Esquivel A., Raga A. C., Cantó J., Rodríguez-González A., López-Cámara D., Velázquez P. F., De Colle F., 2010, ApJ, 725, 1466
  * Gaia Collaboration et al. (2016) Gaia Collaboration et al., 2016, A&A, 595, A1
  * Gaia Collaboration et al. (2018) Gaia Collaboration et al., 2018, A&A, 616, A1
  * Gaia Collaboration et al. (2021) Gaia Collaboration et al., 2021, A&A, 649, A1
  * Garcia-Segura (1997) Garcia-Segura G., 1997, ApJ, 489, L189
  * García-Segura et al. (1999) García-Segura G., Langer N., Różyczka M., Franco J., 1999, ApJ, 517, 767
  * García-Segura et al. (2005) García-Segura G., López J. A., Franco J., 2005, ApJ, 618, 919
  * García-Segura et al. (2016) García-Segura G., Villaver E., Manchado A., Langer N., Yoon S. C., 2016, ApJ, 823, 142
  * Gies & Bolton (1986) Gies D. R., Bolton C. T., 1986, ApJS, 61, 419
  * Gull & Sofia (1979) Gull T. R., Sofia S., 1979, ApJ, 230, 782
  * Gvaramadze & Gualandris (2011) Gvaramadze V. V., Gualandris A., 2011, MNRAS, 410, 304
  * Gvaramadze et al. (2011) Gvaramadze V. V., Röser S., Scholz R. D., Schilbach E., 2011, A&A, 529, A14
  * Gvaramadze et al. (2014) Gvaramadze V. V., Menten K. M., Kniazev A. Y., Langer N., Mackey J., Kraus A., Meyer D. M. A., Kamiński T., 2014, MNRAS, 437, 843
  * Harvey et al. (2020) Harvey E. J., et al., 2020, MNRAS, 499, 2959
  * Henden et al. (2016) Henden A. A., Templeton M., Terrell D., Smith T. C., Levine S., Welch D., 2016, VizieR Online Data Catalog, p. II/336
  * Henney & Arthur (2019) Henney W. J., Arthur S. J., 2019, MNRAS, 486, 3423
  * Herbst et al. (2022) Herbst K., et al., 2022, Space Sci. Rev., 218, 29
  * Hoefner et al. (1996) Hoefner S., Fleischer A. J., Gauger A., Feuchtinger M. U., Dorfi E. A., Winters J. M., Sedlmayr E., 1996, A&A, 314, 204
  * Høg et al. (2000) Høg E., et al., 2000, A&A, 355, L27
  * Huthoff & Kaper (2002) Huthoff F., Kaper L., 2002, A&A, 383, 999
  * Ishihara et al. (2010) Ishihara D., et al., 2010, A&A, 514, A1
  * Kaper et al. (1997) Kaper L., van Loon J. T., Augusteijn T., Goudfrooij P., Patat F., Waters L. B. F. M., Zijlstra A. A., 1997, ApJ, 475, L37
  * Kato (1999) Kato T., 1999, Information Bulletin on Variable Stars, 4789, 1
  * Kobulnicky et al. (2016) Kobulnicky H. A., et al., 2016, ApJS, 227, 18
  * Kobulnicky et al. (2017) Kobulnicky H. A., Schurhammer D. P., Baldwin D. J., Chick W. T., Dixon D. M., Lee D., Povich M. S., 2017, AJ, 154, 201
  * Kukarkin et al. (1970) Kukarkin B. V., Kholopov P. N., Perova N. B., 1970, Information Bulletin on Variable Stars, 480, 1
  * Lamers & Cassinelli (1999) Lamers H. J. G. L. M., Cassinelli J. P., 1999, Introduction to Stellar Winds
  * Langer (2012) Langer N., 2012, ARA&A, 50, 107
  * Lasker et al. (2008) Lasker B. M., et al., 2008, AJ, 136, 735
  * López-Santiago et al. (2012) López-Santiago J., et al., 2012, ApJ, 757, L6
  * McKee (1995) McKee C. F., 1995, in Ferrara A., McKee C. F., Heiles C., Shapiro P. R., eds, Astronomical Society of the Pacific Conference Series Vol. 80, The Physics of the Interstellar Medium and Intergalactic Medium. p. 292
  * Meaburn et al. (2003) Meaburn J., López J. A., Gutiérrez L., Quiróz F., Murillo J. M., Valdéz J., Pedrayez M., 2003, Rev. Mex. Astron. Astrofis., 39, 185
  * Merle et al. (2016) Merle T., Jorissen A., Van Eck S., Masseron T., Van Winckel H., 2016, A&A, 586, A151
  * Meyer et al. (2014) Meyer D. M. A., Mackey J., Langer N., Gvaramadze V. V., Mignone A., Izzard R. G., Kaper L., 2014, MNRAS, 444, 2754
  * Meyer et al. (2016) Meyer D. M. A., van Marle A. J., Kuiper R., Kley W., 2016, MNRAS, 459, 1146
  * Meyer et al. (2017) Meyer D. M. A., Mignone A., Kuiper R., Raga A. C., Kley W., 2017, MNRAS, 464, 3229
  * Meyer et al. (2020) Meyer D. M. A., Oskinova L. M., Pohl M., Petrov M., 2020, MNRAS, 496, 3906
  * Meyer et al. (2021) Meyer D. M. A., Mignone A., Petrov M., Scherer K., Velázquez P. F., Boumis P., 2021, MNRAS, 506, 5170
  * Mohamed et al. (2012) Mohamed S., Mackey J., Langer N., 2012, A&A, 541, A1
  * Neugebauer et al. (1984) Neugebauer G., et al., 1984, ApJ, 278, L1
  * Noriega-Crespo et al. (1997) Noriega-Crespo A., van Buren D., Dgani R., 1997, AJ, 113, 780
  * Olofsson et al. (2010) Olofsson H., Maercker M., Eriksson K., Gustafsson B., Schöier F., 2010, A&A, 515, A27
  * Peri et al. (2012) Peri C. S., Benaglia P., Brookes D. P., Stevens I. R., Isequilla N. L., 2012, A&A, 538, A108
  * Peri et al. (2015) Peri C. S., Benaglia P., Isequilla N. L., 2015, A&A, 578, A45
  * Poveda et al. (1967) Poveda A., Ruiz J., Allen C., 1967, Boletin de los Observatorios Tonantzintla y Tacubaya, 4, 86
  * Raga et al. (1998) Raga A. C., Canto J., Curiel S., Taylor S., 1998, MNRAS, 295, 738
  * Raskin et al. (2011) Raskin G., et al., 2011, A&A, 526, A69
  * Schlegel et al. (1998) Schlegel D. J., Finkbeiner D. P., Davis M., 1998, ApJ, 500, 525
  * Schöier et al. (2006) Schöier F. L., Olofsson H., Lundgren A. A., 2006, A&A, 454, 247
  * Shetye et al. (2018) Shetye S., et al., 2018, A&A, 620, A148
  * Silva Aguirre et al. (2020) Silva Aguirre V., et al., 2020, A&A, 635, A164
  * Skrutskie et al. (2006) Skrutskie M. F., et al., 2006, AJ, 131, 1163
  * Smith (2014) Smith N., 2014, ARA&A, 52, 487
  * Smith & Lambert (1990) Smith V. V., Lambert D. L., 1990, ApJS, 72, 387
  * Spitzer & Varshalovich (1980) Spitzer L., Varshalovich D. A., 1980, Astrofizika, 16, 795
  * Staff et al. (2016) Staff J. E., De Marco O., Macdonald D., Galaviz P., Passy J.-C., Iaconi R., Low M.-M. M., 2016, MNRAS, 455, 3511
  * Steffen et al. (2011) Steffen W., Koning N., Wenger S., Morisset C., Magnor M., 2011, IEEE Transactions on Visualization and Computer Graphics, 17, 454
  * Stevens et al. (1992) Stevens I. R., Blondin J. M., Pollock A. M. T., 1992, ApJ, 386, 265
  * Toalá et al. (2015) Toalá J. A., Guerrero M. A., Ramos-Larios G., Guzmán V., 2015, A&A, 578, A66
  * Villaver & Stanghellini (2005) Villaver E., Stanghellini L., 2005, ApJ, 632, 854
  * Villaver et al. (2012) Villaver E., Manchado A., García-Segura G., 2012, ApJ, 748, 94
  * Vishniac (1994) Vishniac E. T., 1994, ApJ, 428, 186
  * Wareing et al. (2007) Wareing C. J., Zijlstra A. A., O’Brien T. J., Seibert M., 2007, ApJ, 670, L125
  * Weaver et al. (1977) Weaver R., McCray R., Castor J., Shapiro P., Moore R., 1977, ApJ, 218, 377
  * Weigelt et al. (2002) Weigelt G., Balega Y. Y., Blöcker T., Hofmann K. H., Men’shchikov A. B., Winters J. M., 2002, A&A, 392, 131
  * Wenger et al. (2000) Wenger M., et al., 2000, A&AS, 143, 9
  * Wilkin (1996) Wilkin F. P., 1996, ApJ, 459, L31
  * Wilkin & Stahler (1996) Wilkin F. P., Stahler S. W., 1996, in American Astronomical Society Meeting Abstracts. p. 53.05
  * Wright et al. (2010) Wright E. L., et al., 2010, AJ, 140, 1868
  * Zwicky (1957) Zwicky F., 1957, Morphological astronomy
  * van Buren et al. (1995) van Buren D., Noriega-Crespo A., Dgani R., 1995, AJ, 110, 2914
  * van Marle et al. (2011) van Marle A. J., Meliani Z., Keppens R., Decin L., 2011, ApJ, 734, L26
  * van Marle et al. (2012) van Marle A. J., Meliani Z., Marcowith A., 2012, A&A, 541, L8
  * van Marle et al. (2015) van Marle A. J., Meliani Z., Marcowith A., 2015, A&A, 584, A49


11institutetext:  Condensed Matter Physics Laboratory, Heinrich-Heine-
Universität, Universitätsstr.1, 40225 Düsseldorf, Germany

Graphene Electronic transport in mesoscopic systems Scattering by point
defects, dislocations, surfaces, and other imperfections

# Impurity and edge roughness scattering in armchair graphene nanoribbons:
Boltzmann approach

Hengyi Xu  Thomas Heinzel

###### Abstract

The conductivity of armchair graphene nanoribbons in the presence of short-
range impurities and edge roughness is studied theoretically using the
Boltzmann transport equation for quasi-one-dimensional systems. As the number
of occupied subbands increases, the conductivity due to short-range impurities
converges towards the two-dimensional case. Calculations of the
magnetoconductivity confirm the edge-roughness-induced dips at cyclotron radii
close to the ribbon width suggested by the recent quantum simulations.

###### pacs:

81.05.ue

###### pacs:

73.23.-b

###### pacs:

72.10.Fk

##  1 Introduction

The theoretical descriptions of electron transport in graphene sheets with
various scattering sources, such as charged impurities in substrates,
microscopic corrugations or short-range resonant scatters, have been
frequently based on the Boltzmann approach. Such studies, in particular the
electron density dependence of the conductivity, are of fundamental interest
since they help identifying the dominant scattering sources. [1, 2] It has
been verified systematically that the Boltzmann approach works quite well to
describe the transport in broad parameter ranges for both single and double
layer graphene. [3] As the width of the graphene strips is decreased, graphene
nanoribbons (GNRs) are fromed, where size-dependent effects, for example the
inhomogeneous electron density [4] or the edge roughness, become relevant for
the transport properties. In this case, the transport has so far been
described within the framework of the Landauer-Büttiker model with the aid of
Green’s function techniques. [5, 6] It is well established that in GNRs, edge
disorder can contribute significantly to the scattering [7, 8, 9, 10] which in
wide structures is governed by a combination of scattering at charged
impurities and resonant scattering at short-range defects. Edge disorder has
been suggested as the source of the transport gap in narrow GNRs around the
charge neutrality point.[10, 11, 12, 13] Furthermore, a typical size effect,
the so-called edge-roughness-induced magnetoconductance dip (ERID) in GNRs has
been studied by numerical quantum simulations, which is interpreted as a
magnetic-field-enhanced diffusive scattering when the electron trajectory
grazes at the edges. [14] On the other hand, the Boltzmann approach has been
applied to treat a variety of scattering sources in conventional quasi-one-
dimensional (Q1D) systems, for example quantum wires. [15, 16, 17]. However,
only a few aspects of transport in GNRs have so far been studied within the
Boltzmann model. [18, 19]

In the present paper, we apply the linear Boltzmann equation to armchair GNRs
and determine its transport properties in the presence of $\delta$-type short-
range impurities and edge roughness. The magnetoconductivity in wide GNRs with
rough edge roughness is studied as well.

##  2 Model and theory

We start with the Dirac Hamiltonian

$H=\hbar v_{F}(\sigma_{x}\tau_{z}k_{x}+\sigma_{y}k_{y})$ (1)

with Fermi velocity $v_{F}\approx 10^{6}m/s$ and Pauli matrices $\sigma_{x,y}$
and $\tau_{z}$ acting on the $A/B$ sublattice and $K/K^{\prime}$ valley
spaces, respectively. The energy spectrum of GNRs depends on the nature of
their edges, namely zigzag or armchair. Within the present work, we restrict
ourselves to metallic armchair GNRs (AGNRs). For this system, the boundary
conditions imposed on the wave function, namely
$\Psi_{A}(x=0)=\Psi_{B}(x=0)=\Psi_{A}(x=W)=\Psi_{B}(x=W)=0$, give rise to the
allowed transverse wave vectors as

$k_{n}=\frac{n\pi}{W}-\frac{4\pi}{3a}$ (2)

with $a=0.246\mathrm{nm}$ being the lattice constant of graphene. $W$ is the
width of AGNRs. The integer $n$ is of the order of $W/a$ for the energetically
lowest modes. Throughout this text, we denote the energy $\epsilon$ normalized
to $\hbar v_{F}$ as $\widetilde{\epsilon}\equiv\epsilon/(\hbar v_{F})$ with
$\widetilde{\epsilon}^{2}=k_{n}^{2}+k^{2}$. The normalized wave function for
the $n$th subband reads [20, 21]

$\Psi(\mathbf{r})=\frac{e^{iky}}{\sqrt{4WL}}\left(\begin{array}[]{c}e^{ik_{n}x}\\\
\frac{k_{n}+ik}{\widetilde{\epsilon}_{nk}}e^{ik_{n}x}\\\ -e^{-ik_{n}x}\\\ \
-\frac{k_{n}+ik}{\widetilde{\epsilon}_{nk}}e^{-ik_{n}x}\end{array}\right)$ (3)

which is a mixture of two Dirac points $\mathbf{K}=(4\pi/(3a),0)=(K,0)$ and
$\mathbf{K^{\prime}}=(-4\pi/(3a),0)=(-K,0)$. Here we choose the wave vectors
in the $x$-direction to be quantized and the transport is oriented along
$y$-direction. $L$ is the length of the system.

To describe the transport properties of GNRs, we adopt the linearized
Boltzmann equation describing the general Q1D system

$-\frac{eE_{y}}{\hbar}\frac{\partial f_{nk}^{0}(\epsilon_{nk})}{\partial
k}=\sum_{n^{\prime}}\sum_{k^{\prime}}\mathcal{W}_{n^{\prime}k^{\prime}nk}\left[f_{n^{\prime}k^{\prime}}-f_{nk}\right]$
(4)

where $E_{y}$ is the applied electric field along the transport direction,
$f_{nk}$ is the distribution function of a state with wave vector $k$ and
energy $\epsilon_{nk}$ in the $n$th subband, and the superscript “0” denotes
the equilibrium distribution. According to Fermi’s Golden Rule, the scattering
probability due to the perturbation potential is given by

$\mathcal{W}_{n^{\prime}k^{\prime}nk}=\frac{2\pi}{\hbar}|\langle
n^{\prime},k^{\prime}|U|n,k\rangle|^{2}\delta(\epsilon_{n^{\prime}k^{\prime}}-\epsilon_{nk})$
(5)

Using the relaxation time approximation, the nonequilibrium distribution
function can be written as

$f_{nk}(\epsilon_{nk})=f_{nk}^{0}(\epsilon_{nk})-eE_{x}v_{n}(\epsilon_{nk})\tau_{n}(\epsilon_{nk})\delta(\epsilon_{nk}-E_{F})$
(6)

with Fermi energy $E_{F}$ and the relaxation time $\tau_{n}$ for the state in
the $n$th subband. The velocity for the $n$th subband
$v_{n}=(1/\hbar){\partial\epsilon_{nk}}/{\partial
k}=v_{F}{k}/{\sqrt{(k_{n})^{2}+k^{2}}}$ for the linear spectrum of graphene.
Inserting Eq. (6) into Eq. (4), the Boltzmann equation at zero temperature can
be written as

$\displaystyle\frac{k}{\widetilde{\epsilon}_{nk}}\delta(\epsilon_{nk}-E_{F})$
$\displaystyle=$
$\displaystyle\sum_{n^{\prime},k^{\prime}}\mathcal{W}_{n^{\prime}k^{\prime}nk}\left[\frac{k}{\widetilde{\epsilon}_{nk}}\tau_{n}(\epsilon_{nk})\delta(\epsilon_{nk}-E_{F})\right.$
(7)
$\displaystyle\left.-\frac{k^{\prime}}{\widetilde{\epsilon}_{n^{\prime}k^{\prime}}}\tau_{n^{\prime}}(\epsilon_{n^{\prime}k^{\prime}})\delta(\epsilon_{n^{\prime}k^{\prime}}-E_{F})\right].$

Multiplying both sides of Eq. (7) by $k$ and summing over $k$, we obtain after
some algebra

$k_{F}^{n}=\frac{2\pi}{\hbar}\sum_{n^{\prime}}\mathcal{T}_{nn^{\prime}}\tau_{n^{\prime}}(E_{F})$
(8)

where $k_{F}^{n}$ is Fermi wave vector in the $n$th subband. The transition
matrix element $\mathcal{T}_{nn^{\prime}}$ is defined as

$\displaystyle\mathcal{T}_{nn^{\prime}}$ $\displaystyle=$
$\displaystyle\frac{\pi\hbar
v_{F}}{L}\sum_{k^{\prime}}\sum_{k}\left[\delta_{nn^{\prime}}\sum_{\mu}|\langle\mu,k^{\prime}|U|n,k\rangle|^{2}\right.$
$\displaystyle\times$
$\displaystyle\frac{k^{2}}{\widetilde{\epsilon}_{nk}}\delta(\epsilon_{nk}-E_{F})\delta(E_{\mu
k^{\prime}}-E_{F})$ $\displaystyle-$ $\displaystyle\left.|\langle
n^{\prime},k^{\prime}|U|n,k\rangle|^{2}\frac{kk^{\prime}}{\widetilde{\epsilon}_{n^{\prime}k^{\prime}}}\delta(\epsilon_{nk}-E_{F})\delta(\epsilon_{n^{\prime}k^{\prime}}-E_{F})\right]$

with the summation running over the mode index $\mu$. The Boltzmann
conductivity for GNRs is then given by

$\sigma(E_{F})=\frac{2e^{2}}{h}\frac{\hbar^{2}v_{F}^{2}}{\pi
E_{F}}\frac{1}{W}\sum_{n,n^{\prime}}k_{F}^{n}k_{F}^{n^{\prime}}(\mathcal{T}^{-1})_{nn^{\prime}}$
(10)

For nonzero temperature, the conductivity is obtained from

$\sigma=\int d\epsilon\left(-\frac{\partial
f(\epsilon)}{\partial\epsilon}\right)\sigma(\epsilon)$ (11)

We proceed by describing how the scattering potentials have been implemented
in this formalism. First, we consider $\delta$-type impurities in the form of

$U=\gamma\sum_{j=1}^{N_{I}}\delta(x-x_{j})\delta(y-y_{j})$ (12)

where $\gamma$ and $N_{I}$ are the strength and the number of impurities,
respectively. Thus the matrix element squared of the perturbation is evaluated
as

$\displaystyle|\langle
n^{\prime}k^{\prime}|U|nk\rangle|^{2}=\frac{\gamma^{2}}{4W^{2}L^{2}}$ (13)
$\displaystyle\times$
$\displaystyle\sum_{j=1}^{N_{j}}\cos^{2}\frac{(n-n^{\prime})\pi
x_{j}}{W}\left|1+\frac{(k_{n^{\prime}}-ik^{\prime})(k_{n}+ik)}{\widetilde{\epsilon}_{n^{\prime}k^{\prime}}\widetilde{\epsilon}_{nk}}\right|^{2}$
$\displaystyle=$
$\displaystyle\frac{\gamma^{2}n_{i}}{4WL}(1+\delta_{nn^{\prime}})\left(1+\frac{k_{n^{\prime}}k_{n}+k^{\prime}k}{\widetilde{\epsilon}_{n^{\prime}k^{\prime}}\widetilde{\epsilon}_{nk}}\right)$

with $n_{i}=N_{I}/WL$. The transition matrix elements are finally given by

$\displaystyle\mathcal{T}_{nn^{\prime}}$ $\displaystyle=$
$\displaystyle\frac{\gamma^{2}n_{i}}{4\pi
W}\delta_{nn^{\prime}}\left(\sum_{\mu}\left[(1+\delta_{n\mu})\left(\frac{E_{F}}{\hbar^{2}v_{F}^{2}}+\frac{k_{\mu}k_{n}}{E_{F}}\right)\frac{k_{F}^{n}}{k_{F}^{\mu}}\right]\right.$
(14)
$\displaystyle\left.-\frac{k_{F}^{n^{\prime}}k_{F}^{n}}{{E}_{F}}\right)-\frac{\gamma^{2}n_{i}}{4\pi
W}\frac{k_{F}^{n^{\prime}}k_{F}^{n}}{{E}_{F}}.$

Using Eq. (14) and Eq. (10) we obtain for the conductivity of GNRs the
expression

$\displaystyle\sigma$ $\displaystyle=$
$\displaystyle\frac{8e^{2}}{h}\frac{(\hbar
v_{F})^{2}}{\gamma^{2}n_{i}}\sum_{n,n^{\prime}}\left\\{\left[\sum_{\mu}\left((1+\delta_{n\mu})(\widetilde{E}^{2}_{F}+k_{\mu}k_{n})\frac{k_{F}^{n}}{k_{F}^{\mu}}\right)\right.\right.$
(15)
$\displaystyle\left.\left.-k_{F}^{n^{\prime}}k_{F}^{n}\right]\delta_{nn^{\prime}}-k_{F}^{n^{\prime}}k_{F}^{n}\right\\}^{-1}_{nn^{\prime}}k_{F}^{n}k_{F}^{n^{\prime}}$

For a large number of subbands, i.e., $\mathcal{N}\gg 1$, we have
$k_{\mu},k_{n}\ll\widetilde{E}_{F}$ and
$k_{F}^{n},k_{F}^{n^{\prime}}\sim\widetilde{E}_{F}$. In this approximation,
Eq. (15) converges to the well-known result for the extended case

$\sigma=\frac{8e^{2}}{h}\frac{\hbar^{2}v_{F}^{2}}{\gamma^{2}n_{i}}$ (16)

which is independent of the carrier concentration. [2, 3]

As a second scattering mechanism, we study the effects of edge roughness in
the absence of magnetic fields. The edge roughness is parametrized by
$\Delta(y){\partial V(x)}/{\partial x}$, an expression which has been applied
before to model rough semiconductor quantum wires and interfaces. [22, 23]
$V(x)$ is the one-dimensional confinement potential which can be modeled by a
finite mass term in the Dirac Hamiltonian. [24] $\Delta(y)$ is a function
describing the potential fluctuation of GNRs and characterized by
$\langle\Delta(y)\rangle=0$ and the autocovariance function
$\langle\Delta(y)\Delta(y^{\prime})\rangle=\Delta^{2}\exp[-(y-y^{\prime})^{2}/\Lambda^{2}]$
with $\Lambda$ being the correlation length. Furthermore,
$\langle\cdots\rangle$ denotes position averaging.

To evaluate the perturbation matrix element, we define the function $\Xi$
related to the $x$-components of the wave functions as

$\Xi_{n^{\prime}n}=\frac{1}{W}\int_{0(W)}^{-\infty(+\infty)}dx\phi_{n^{\prime}}^{*}(x)\frac{\partial
V(x)}{\partial x}\phi_{n}(x)$ (17)

where $\phi_{n}(x)$ denotes one of the components of wave function in Eq. (3).
For the hard-wall confinement potential present in GNRs, this function can be
expressed as

$\Xi_{n^{\prime}n}=\left.-\frac{1}{W}\frac{\hbar
v_{F}}{2\widetilde{E}_{F}}\left[\frac{\partial\phi^{*}_{n^{\prime}}}{\partial
x}\frac{\partial\phi_{n}}{\partial x}\right]\right|_{x=0,W}$ (18)

It is noteworthy that a linear form for the matrix elements of the edge
roughness perturbation has been used, which however neglects the interband
scattering. [25, 19] Using Eq. (18), the square of the matrix element for edge
roughness reads

$\displaystyle|\langle n^{\prime},k^{\prime}|U|n,k\rangle|^{2}$
$\displaystyle=$ $\displaystyle\frac{\pi^{9/2}n^{\prime
2}n^{2}}{8W^{6}}\frac{(\hbar
v_{F})^{2}}{\widetilde{E}^{2}_{F}}\left(1+\frac{k_{n}k_{n^{\prime}}+kk^{\prime}}{\widetilde{\epsilon}_{nk}\widetilde{\epsilon}_{n^{\prime}k^{\prime}}}\right)$
(19)
$\displaystyle\times\frac{\Lambda\Delta^{2}}{L}\exp[-\Lambda^{2}(k-k^{\prime})^{2}/4]$

where we have used the Gaussian integral in the evaluation of part in the
$y$-direction. In the case of small correlation length
$\Lambda\ll\lambda_{F}$, the transition matrix element has the form

$\displaystyle\mathcal{T}_{nn^{\prime}}$ $\displaystyle=$
$\displaystyle\frac{\pi^{7/2}\Lambda\Delta^{2}}{8W^{6}}\frac{\hbar
v_{F}}{\widetilde{E}_{F}}n^{2}\left[\sum_{\mu}\mu^{2}\left(1+\frac{k_{n}k_{\mu}}{\widetilde{E}^{2}_{F}}\right)\frac{k_{F}^{n}}{k_{F}^{\mu}}\delta_{nn^{\prime}}\right.$
(20) $\displaystyle\left.-n^{\prime
2}\frac{k_{F}^{n^{\prime}}k_{F}^{n}}{\widetilde{E}^{2}_{F}}\right]$

This results in a conductivity given by

$\displaystyle\sigma$ $\displaystyle=$
$\displaystyle\frac{16e^{2}}{h}\frac{W^{5}}{\pi^{9/2}\Lambda\Delta^{2}}\sum_{n,n^{\prime}}k_{F}^{n}k_{F}^{n^{\prime}}$
$\displaystyle\times$
$\displaystyle\left[n^{2}\sum_{\mu}\mu^{2}\left(1+\frac{k_{n}k_{\mu}}{\widetilde{E}^{2}_{F}}\right)\frac{k_{F}^{n}}{k_{F}^{\mu}}\delta_{nn^{\prime}}-n^{2}n^{\prime
2}\frac{k_{F}^{n^{\prime}}k_{F}^{n}}{\widetilde{E}^{2}_{F}}\right]^{-1}_{nn^{\prime}}.$

Regarding the diagonal contributions of the transition rate matrix, it is
convenient to write the inverse of the relaxation time for the $n$th subband
as

$\frac{1}{\tau_{n}}=\frac{\pi^{9/2}}{4W^{6}}\frac{\hbar
v_{F}^{2}}{E_{F}}\Lambda\Delta^{2}n^{2}\sum_{\mu}\mu^{2}\left[\frac{1}{k_{F}^{\mu}}+\frac{k_{\mu}k_{n}/k_{F}^{\mu}}{\widetilde{E}^{2}_{F}}\right].$
(22)

For a large number of occupied subbands ($\mathcal{N}\gg 1$), the second term
in the bracket can be neglected since $k_{n,\mu}\ll\widetilde{E}_{F}$, and the
summation can be replaced by an integral. Eq.(22) can then be written as

$\frac{1}{\tau_{n}}\approx\frac{\pi^{5/2}}{16W^{3}}\Lambda\Delta^{2}v_{F}\widetilde{E}_{F}n^{2}$
(23)

which shows a striking similarity for the corresponding results reported for
semiconductor quantum wires. [15]. Furthermore, this relaxation time results
in a conductivity in the limit of $\mathcal{N}\gg 1$ of

$\sigma\approx\frac{32e^{2}}{h}\frac{1}{3\pi^{1/2}}\frac{W^{2}}{\Lambda\Delta^{2}\widetilde{E}_{F}}.$
(24)

Up to now, we have looked at the transport with edge roughness in the absence
of magnetic fields $B$ and now continue by including it and discussing its
semiclassical effects. A weak magnetic field may homogenize the contributions
of the occupied subbands to the overall conductivity and result in a reduction
of the magnetoconductivity when the cyclotron radius $r_{c}\sim W$. This
effect has been verified numerically by recent quantum calculations in GNRs
[14] but, to the best of our knowledge, has so far not yet been observed. Here
the maximum reduction of the relaxation time $\tau$ in magnetic fields can be
roughly estimated by averaging over all occupied modes, resulting in

$\frac{1}{\tau(B>0)}\approx\frac{\pi^{1/2}}{48W}\Lambda\Delta^{2}v_{F}\widetilde{E}^{3}_{F}.$
(25)

The magnetoconductivity near the dip is then given by

$\sigma(B>0)\approx\frac{48e^{2}}{h}\frac{1}{\pi^{1/2}}\frac{W}{\Lambda\Delta^{2}\widetilde{E}_{F}^{2}},$
(26)

which holds for $\mathcal{N}\gg 1$.

##  3 Numerical results and discussion

Figure 1: (Colour online) The numerical conductivity of $\delta$-type
impurities in armchair GNRs with $W=180$nm plotted versus the Fermi energy for
two different temperatures T. We have chosen the parameters as
$\overline{\gamma}=0.9$ and $\overline{n}_{i}=0.5$, where
$\overline{\gamma}=\gamma(\hbar v_{F})^{2}/E_{F}$ and
$\overline{n}_{i}=n_{i}\lambda_{F}^{2}$.

In Fig. 1 we show the conductivity for $\delta$-type impurities according to
Eq.(15) as a function of Fermi energy. For different degrees of disorder, the
conductivity shows very similar features while the amplitude of the
conductivity depends on the disorder parameters. Prominent quantum
oscillations at zero temperature are observed, i.e., the conductivity drops
rapidly as a new scattering channel is opened and increases again until the
Fermi energy hits the next subband at larger energies. The oscillations are
smeared by finite temperature to some extent. In the whole range of Fermi
energies studied, the average conductivity remains independent of the carrier
concentration, which is consistent with the two-dimensional case. [2, 3]

Figure 2: (Colour online) The Boltzmann conductivity as a function of the
Fermi energy for different edge roughness in AGNRs. The dashed and solid lines
correspond to zero temperature and a temperature of $T=10\mathrm{K}$,
respectively. The smooth solid lines are calculated from Eq.(24) in the limit
of $\mathcal{N}\gg 1$.

Fig. 2 shows the conductivity for edge roughness as a function of Fermi energy
calculated from Eq. (LABEL:sigma_edge). The parameter values chosen for
$\Lambda$ and $\Delta$ correspond to short-range defects, for instance, a few
atoms missing at the GNR edges, as widely assumed in simulations of edge
disorder [8, 26, 9, 14]. The correlation length ensures that
$\Lambda\ll\lambda_{F}$ over the whole range of Fermi energies. The Boltzmann
conductivity at nonzero temperature (indicated by the solid lines) shows
suppressed quantum fluctuations in comparison with the zero-temperature cases
(indicated by the dashed line). The overall conductivity decreases as the
Fermi energy increases. Since the correlation length $\Lambda$ and edge
position fluctuation amplitude $\Delta$ increase relative to the Fermi
wavelength as $E_{F}$ is increased, this behavior is similar to that one found
in conventional quantum wires [15, 17]. In the case of large number of
subbands, $\mathcal{N}\gg 1$, the results from Eq. (24) (indicated by solid
lines) exhibit the same overall trends and agree well with the exact ones
except the absence of the quantum oscillations.

The Fermi energies in Fig. 1 and 2 correspond to numbers of subbands between
$10$ and $30$. For smaller Fermi energies, i.e. a few occupied modes,
conductivity shows more prominent quantum fluctuations and may deviate
considerably from the asymptotic expressions. Moreover, it should be noted
that our calculations based on the Boltzmann approach is valid for the case
where the interband scattering is rather strong. This is the case when the
mean free path is considerably shorter than the length of the graphene ribbon.

Figure 3: (Colour online) The width dependence of the conductivity for
different edge roughness in AGNRs with $E_{F}=200\mathrm{meV}$ and temperature
$T=10\mathrm{K}$. The solid lines show the exact conductivity from Eq.(11) by
the $\mathcal{T}$ matrix inversion, and dashed lines correspond to the limit
$\mathcal{N}\gg 1$ from Eq.(24).

Fig. 3 shows the Boltzmann conductivity as a function of the GNR width for
different edge roughness parameters. Here, only results for $T>0$ are
presented. The overall conductivity for two roughness levels exhibits a
parabolic dependence on the width, superimposed by quantum oscillations. This
quadratic behavior may be seen more clearly from the analytical expression Eq.
(24), as illustrated by the dashed lines.

Figure 4: (Colour online) The magnetoconductivity around ERID as a function of
the width for different Fermi energies at finite temperature $T=10\mathrm{K}$.
(Note the logarithmic scale.) The zero-field conductivity with
$E_{F}=200\mathrm{meV}$ is also plotted at the top. The solid and dashed lines
correspond to the results from Eq.(LABEL:sigma_edge) and Eq.(26) for
$\mathcal{N}\gg 1$, respectively. $r_{c}$ is the cyclotron radius.

In the following, we give a rough estimate of the GNR conductivity in magnetic
fields with amplitudes close to the position of the ERID, i.e., $r_{c}\approx
W$. A more exact calculation would have to rely upon a calculation of the wave
functions in magnetic fields, which can be obtained by solving the
eigenequation of the Dirac Hamiltonian with magnetic fields included. [27, 28]
This, however, should have only a marginal effect and we limit ourselves to
the qualitative properties of the system close to the ERID.

In Fig. 4, the conductivity at $T=10\mathrm{K}$ is shown in the vicinity of
the ERID for different Fermi energies. The parameters for edge roughness are
fixed to $\Lambda=0.3\mathrm{nm}$ and $\Delta=6\mathrm{nm}$. For comparison,
the corresponding zero-field conductivity with $E_{F}=200\mathrm{meV}$ is
plotted as well. The conductivity around the ERID increases linearly with the
GNR width, in contrast to the parabolic dependence in the absence of a
magnetic field. This linear relationship can be easily identified from Eq.
(26) and is also illustrated by the dashed lines in Fig. 4. As a consequence
of the distinctly different dependencies of $\sigma$ on W, the ERID can be
expected to be more pronounced in wider GNRs. This feature is in qualitative
agreement with our previous quantum simulations. [14].

We conclude this analysis by commenting on the observability of the ERID in
realistic GNRs. First of all, the length of the GNR is irrelevant in the
present treatment since diffusive transport has been assumed. Second, we have
restricted ourselves to the case of rather small correlation lengths for the
edge roughness. It is self-evident that a large correlation length suppresses
the ERID in view the reduced diffusiveness of the scattering at the edges.
Furthermore, the bulk disorder must remain at a sufficiently low level as
indicated by the quantum simulations before, such that it does not mask
completely the edge roughness scattering. We moreover expect that
qualitatively, the ERID does not depend much in the the type of edges, even
though numerical simulations suggest that zigzag GNRs are more robust with
respect to edge disorder. [14] Similar analytical expressions for zigzag GNRs
are possible in principle but more complicated due to the presence of surface
states and the interdependence of the transverse and longitudinal wave
vectors.

In summary, we have studied the transport properties of AGNRs with short-range
impurities and edge roughness within the framework given by the Boltzmann
equation. An edge-roughness-induced magnetoconductivity minimum suggested by
the recent quantum calculations is confirmed by the Boltzmann results and
should become observable experimentally if the correlation length of the edge
roughness is not much larger than the Fermi wavelength and the bulk disorder
is sufficiently low. It has been shown that the ERID induced by the magnetic-
field-enhanced diffusive scattering at rough edges shows a behavior very
similar to that one found in conventional semiconductor quantum wires, despite
the fundamentally different energy dispersion.

###### Acknowledgements.

H.X. and T.H. acknowledge financial support from Heinrich-Heine-Universität
Düsseldorf.

## References

  * [1] Peres N. M. R. Rev. Mod. Phys.8220112673.
  * [2] Sarma S. D., Adam S., Hwang E. H. Rosso E. Rev. Mod. Phys.832011407.
  * [3] Xu H., Heinzel T. Zozoulenko I. V. Phys. Rev. B842011115409.
  * [4] Xu H., Heinzel T., Shylau A. A. Zozoulenko I. V. Phys. Rev. B822010115311.
  * [5] Lewenkopf C. H., Mucciolo E. R. Castro Neto A. H. Phys. Rev. B772008081410.
  * [6] Xu H., Heinzel T., Evaldsson M. Zozoulenko I. V. Phys. Rev. B772008245401.
  * [7] Areshkin D. A., Gunlycke D. White C. T. Nano Lett.72007204.
  * [8] Evaldsson M., Ihnatsenka S. Zozoulenko I. V. Phys. Rev. B772008165306.
  * [9] Xu H., Heinzel T. Zozoulenko I. V. Phys. Rev. B802009045308.
  * [10] Han M. Y., Ozyilmaz B., Zhang Y. Kim P. Phys. Rev. Lett.982007206805.
  * [11] Stampfer C., Guettinger J., Hellmueller S., Molitor F., Ensslin K. Ihn T. Phys. Rev. Lett.1022009056403.
  * [12] Liu X. L., Oostinga J. B., Morpurgo A. F. Vandersypen L. M. K. Phys. Rev. B802009121407.
  * [13] Todd K., Chou H. T., Amasha S. Goldhaber-Gordon D. Nano Lett.92009416.
  * [14] Xu H., Heinzel T. Zozoulenko I. V. EPL (Europhysics Letters)97201228008.
  * [15] Akera H. Ando T. Phys. Rev. B11676199143.
  * [16] Bruus H., Flensberg K. Smith H. Phys. Rev. B48199311144.
  * [17] Feilhauer J. Moško M. Phys. Rev. B832011245328.
  * [18] Huang D. Gumbs G. J. Appl. Phys.1072010103710.
  * [19] Huang D., Gumbs G. Roslyak O. Phys. Rev. B832011115405.
  * [20] Neto A. H. C., Guinea F., Peres N. M. R., Novoselov K. S. Geim A. K. Rev. Mod. Phys.812009109.
  * [21] Wurm J., Wimmer M., Adagideli ., Richter K. Baranger H. U. New Journal of Physics112009095022.
  * [22] Goodnick S. M., Ferry D. K., Wilmsen C. W., Liliental Z., Fathy D. Krivanek O. L. Phys. Rev. B3219858171.
  * [23] Ferry D. K. Goodnick S. M. Transport in Nanostructures 1st Edition (Cambridge University Press) 1997.
  * [24] Peres N. M. R., Rodrigues J. N. B., Stauber T. dos Santos J. M. B. L. Journal of Physics: Condensed Matter212009344202.
  * [25] Fang T., Konar A., Xing H. Jena D. Phys. Rev. B782008205403.
  * [26] Mucciolo E. R., Castro Neto A. H. Lewenkopf C. H. Phys. Rev. B792009075407.
  * [27] Brey L. Fertig H. A. Phys. Rev. B732006195408.
  * [28] De Martino A., Hütten A. Egger R. Phys. Rev. B842011155420.