# Synthetic observations of molecular clouds in a galactic center environment:
I. Studying maps of column density and integrated intensity

Erik Bertram1, Simon C. O. Glover1, Paul C. Clark2, Sarah E. Ragan3 & Ralf S.
Klessen1
1Universität Heidelberg, Zentrum für Astronomie, Institut für Theoretische
Astrophysik, Albert-Ueberle-Str. 2, 69120 Heidelberg, Germany
2School of Physics and Astronomy, Cardiff University, CF24 3AA, UK
3School of Physics and Astronomy, University of Leeds, LS2 9JT, England

###### Abstract

We run numerical simulations of molecular clouds (MCs), adopting properties
similar to those found in the Central Molecular Zone (CMZ) of the Milky Way.
For this, we employ the moving mesh code Arepo and perform simulations which
account for a simplified treatment of time-dependent chemistry and the non-
isothermal nature of gas and dust. We perform simulations using an initial
density of $n_{0}=10^{3}\,$cm-3 and a mass of $1.3\times 10^{5}\,$M⊙.
Furthermore, we vary the virial parameter, defined as the ratio of kinetic and
potential energy, $\alpha=E_{\text{kin}}/|E_{\text{pot}}|$, by adjusting the
velocity dispersion. We set it to $\alpha=0.5,2.0$ and $8.0$, in order to
analyze the impact of the kinetic energy on our results. We account for the
extreme conditions in the CMZ and increase both the interstellar radiation
field (ISRF) and the cosmic-ray flux (CRF) by a factor of 1000 compared to the
values found in the solar neighbourhood. We use the radiative transfer code
Radmc-3d to compute synthetic images in various diagnostic lines. These are
[Cii] at 158 $\mu$m, [Oi] (145 $\mu$m), [Oi] (63 $\mu$m), 12CO
($J=1\rightarrow 0$) and 13CO ($J=1\rightarrow 0$) at 2600 $\mu$m and 2720
$\mu$m, respectively. When $\alpha$ is large, the turbulence disperses much of
the gas in the cloud, reducing its mean density and allowing the ISRF to
penetrate more deeply into the cloud’s interior. This significantly alters the
chemical composition of the cloud, leading to the dissociation of a
significant amount of the molecular gas. On the other hand, when $\alpha$ is
small, the cloud remains compact, allowing more of the molecular gas to
survive. We show that in each case the atomic tracers accurately reflect most
of the physical properties of both the H2 and the total gas of the cloud and
that they provide a useful alternative to molecular lines when studying the
ISM in the CMZ.

###### keywords:

Galaxy: centre – galaxies: ISM – ISM: clouds – stars: formation

##  1 Introduction

The inner few hundred parsecs of our Milky Way, known as the Central Molecular
Zone (CMZ), are rich in molecular and dense gas and account for about 5% of
the molecular gas content of the Galaxy (Morris & Serabyn, 1996). As shown by
Longmore et al. (2013), the CMZ contains gas with densities of the order of
several $10^{3}\,$cm-3 and a total mass of $\sim 10^{7}\,$M⊙. Furthermore, it
is highly turbulent (Shetty et al., 2012) and illuminated by a strong
interstellar radiation field (ISRF) combined with a high cosmic-ray flux (CRF)
(Yusef-Zadeh et al., 2007; Clark et al., 2013). It thus provides an excellent
laboratory for studying the physics of the interstellar medium (ISM) under
extreme conditions (such as those that also occur in distant starburst
galaxies) with very high resolution.

Previous observations suggested the existence of copious CO gas in the CMZ
near the Galactic Center (GC) (Bania, 1977; Burton et al., 1978; Liszt &
Burton, 1978; Morris & Serabyn, 1996; Bitran et al., 1997; Oka et al., 1998;
Martin et al., 2004), which might be used to trace the density and velocity
structure of the underlying molecular component of the ISM. However, the
effect of the extreme physical conditions in the CMZ on the CO distribution
within the gas and the degree to which CO is a biased tracer of the underlying
cloud properties remain relatively unexplored issues. Previous numerical
models of clouds in the CMZ (Clark et al., 2013) or in other harsh
environments (Clark & Glover, 2015) have shown that the gas is likely to be
highly chemically inhomogeneous. This means that molecular tracers might tell
us only little about the kinematics, the temperature and density distributions
of the cloud and so on. Hence, in order to learn more about the internal
physics (e.g. heating, cooling and the chemistry) of the warm gas, which makes
up a large fraction of the total mass, we have to use other probes.

Recently, Clark et al. (2013) have shown that the [Oi] 63 $\mu$m and the [Cii]
158 $\mu$m fine structure lines dominate the cooling in CMZ clouds over a wide
range in densities. Moreover, Rodríguez-Fernández et al. (2004) present
observations of these fine structure lines in the CMZ, carried out with the
ISO satellite. They show that both lines are very strong in this region and
can be used to infer important information about the physics of the gas in the
CMZ. Hence, using such atomic tracers in order to study the physics of MCs in
the CMZ seems to be a promising approach.

In this paper we follow up on this idea and investigate the applicability of
atomic tracers to study the properties of the ISM in the CMZ. We aim to
provide synthetic observations, readily comparable to observations of fine
structure lines. To do so, we compute synthetic observations of MCs in a CMZ-
like environment using the radiative transfer code Radmc-3d (Dullemond, 2012)
for various diagnostic lines, specifically for [Cii] (158 $\mu$m), [Oi] (145
$\mu$m), [Oi] (63 $\mu$m), 12CO (2600 $\mu$m) and 13CO (2720 $\mu$m). We
perform simulations of various clouds with the moving mesh code Arepo
(Springel, 2010) using environmental properties similar to those experienced
by a typical CMZ cloud. In particular, we adopt values for the ISRF strength
and the CRF comparable to those inferred for the dense CMZ cloud G0.253+0.016,
also known as “The Brick” (see e.g. Güsten et al., 1981; Lis et al., 1994; Lis
& Menten, 1998; Lis et al., 2001; Molinari et al., 2011; Immer et al., 2012;
Longmore et al., 2012; Kauffmann et al., 2013; Clark et al., 2013; Johnston et
al., 2014; Rathborne et al., 2014, 2015). We use values for the ISRF and the
CRF that mimic the harsh conditions assumed to be found in the CMZ, which are
a factor of $\sim 1000$ larger than the values measured in the solar
neighbourhood (Habing, 1968; Draine, 1978; Mathis et al., 1983). We simulated
clouds with an initial number density of $n_{0}=10^{3}\,$cm-3 and studied the
impact of different virial $\alpha$ parameters, $\alpha=0.5,2.0$ and $8.0$, on
the physical properties. More information about the simulations can be found
in Bertram et al. (2015). Furthermore, we note that this study is the first of
two papers which analyze the properties of CMZ-like clouds. In this paper, we
focus on maps of column density and integrated intensity, while in Paper II
(Bertram et al., in prep.), we study the kinematic properties of our model
clouds.

This paper is structured as follows. In section 2 we present our numerical
simulations and the radiative transfer post-processing tool. In section 3 we
show and discuss the results of our studies. We present a summary and our
conclusions in section 4.

##  2 Methods

Model name | Virial $\alpha$ | Initial number density $n_{0}$ | 3D velocity dispersion $\sigma_{v}$ | Radius $R$ | Free-fall time $t_{\text{ff}}$ | Simulation end $t_{\text{end}}$
---|---|---|---|---|---|---
|  | [cm-3] | [km/s] | [pc] | [Myr] | [Myr]
GC-0.5-1000 | 0.5 | 1000 | 5.4 | 8.9 | 1.40 | 1.00
GC-2.0-1000 | 2.0 | 1000 | 10.8 | 8.9 | 1.40 | 1.40
GC-8.0-1000 | 8.0 | 1000 | 21.6 | 8.9 | 1.40 | 1.40
Table 1: Overview of the initial conditions for our different cloud models.

The numerical simulations studied in this paper are described in more detail
in Bertram et al. (2015). However, we summarize the most important aspects of
the simulations here and introduce the radiative transfer post-processing,
which we use for computing synthetic observational maps of the clouds in
various diagnostic lines.

###  2.1 Hydrodynamical and chemical model

We use the moving mesh code Arepo (Springel, 2010) in order to run numerical
simulations of various model clouds. We have added to Arepo a detailed atomic
and molecular cooling function, described in detail in Glover et al. (2010)
and Glover & Clark (2012), and a simplified treatment of the basic chemistry
of the gas. The chemical network is based on the work of Nelson & Langer
(1997) and Glover & Mac Low (2007), and allows us to follow the formation and
destruction of H2 and CO self-consistently within our simulations. The network
tracks the abundances of 6 species and follows 14 chemical reactions. The
simplified Nelson & Langer (1997) network is known to somewhat overestimate
the rate at which CO forms, owing to its neglect of atomic carbon (see the
detailed discussion in Glover & Clark, 2012). It also neglects the reaction

$\text{CO}+\text{He}^{+}\rightarrow\text{C}^{+}+\text{O}+\text{He},$ (1)

which plays an important role in regulating the CO abundance in gas exposed to
a high CRF (see, e.g. Clark & Glover, 2015; Bisbas et al., 2015). Our
simulations thereby overestimate the CO abundance in the simulated clouds, and
hence potentially underestimate the C+ and O abundance in regions where the CO
abundance is large. Therefore, in order to establish the effect that this has
on the synthetic maps, we make two versions of each map: one in which we use
values for the fractional abundances of C+ and O taken from the simulations
and a second version in which we assume that all of the carbon is in the form
of C+ and all of the oxygen is present as O. This can be achieved by assuming
zero CO abundances in our simulations. As discussed in more detail in Appendix
C, we find only minor differences in the [Oi] (145 $\mu$m) emission PDF and
negligible differences in the [Oi] (63 $\mu$m) and [Cii] (158 $\mu$m) PDFs. We
therefore conclude that the known weaknesses of our simplified chemical model
have little influence on the results we obtain for [Cii] and [Oi]. Full
details of the chemical model with a description of how the chemistry
interacts with the ISRF via the TreeCol algorithm can be found in Clark et al.
(2012). Further examples of the use of our chemical model with the Arepo code
can be found in Smith et al. (2014a) and Smith et al. (2014b).

We assume that the gas has a uniform solar metallicity and adopt the standard
ratio of helium to hydrogen. The abundances of carbon and oxygen are taken
from Sembach et al. (2000). We use $x_{\text{C}}=1.4\times 10^{-4}$ and
$x_{\text{O}}=3.2\times 10^{-4}$, where $x_{\text{C}}$ and $x_{\text{O}}$ are
the fractional abundances by number of carbon and oxygen relative to hydrogen.
We note that the CMZ has super-solar metallicity. Nevertheless, we use a
uniform solar value in order to be conservative regarding the cooling rates in
our runs. When we start the simulations, hydrogen, helium and oxygen are in
atomic form, while carbon is assumed to be in singly ionized form, as C+. We
adopt the standard local value for the dust-to-gas ratio of 1:100 (for further
discussion see, e.g. Glover et al., 2010). We set the cosmic ray ionization
rate (CRIR) of atomic hydrogen to $\zeta=3\times 10^{-14}\>{\rm s^{-1}}$
(Clark et al., 2013). This value is a factor $\sim 1000$ higher than the value
measured in dense clouds in the solar neighbourhood (van der Tak & van
Dishoeck, 2000) and is comparable to the high value that has previously been
inferred in the CMZ (Yusef-Zadeh et al., 2007; Clark et al., 2013). Note,
however, that it is only around a factor of 100 higher than the value derived
in more diffuse gas in the solar neighbourhood (Indriolo et al., 2015). For
the interstellar radiation field, we adopt the same spectral shape as given in
Draine (1978). We denote the strength of the Draine ISRF as $G_{0}=1$ and
perform simulations with a field strength $G_{0}=1000$ (Clark et al., 2013).
This corresponds to an integral flux in the energy range $6-13.6\,$eV of
$2.7\times 10^{-3}\,$erg cm-2s-1.

###  2.2 Model parameters

The clouds that we model are initially spherical and located at the center of
a large box of low-density gas. For the box, we use periodic boundary
conditions. However, this is simply a convenient choice, since the boundaries
do not influence the dense cloud evolution, because the size of the box is
taken to be much larger than the size of the cloud. The cloud has a uniform
initial hydrogen nuclei number density, which we set to $10^{3}\,$cm-3, and a
total mass of $M_{\text{tot}}=1.3\times 10^{5}\,$M⊙. The initial cloud radius
is $R\approx 8.9\,$pc. In all of our simulations, the density of the gas
surrounding the cloud is $\approx 1$ cm-3, but we note that our results are
insensitive to this value provided that it is much smaller than the mean cloud
density. The cubic side length of the total simulation domain is set to
$5\times$ the individual cloud radius of 8.9 pc, i.e. 44.5 pc in total. All
clouds have zero bulk velocity and are placed in the center of the box. The
Voronoi cells initially have approximately constant volumes. We initially
start with $2\times 10^{6}$ cells in total. The initial cell mass within the
cloud corresponds to $\approx 2\,$M⊙.

We make use of a Jeans refinement criterion in order to accurately refine
dense and collapsed gas regions in the box over the whole simulation time. We
use a constant number of 8 cells per Jeans length, which is sufficient to
avoid artificial fragmentation (see, e.g. Truelove et al. 1998, Greif et al.
2011 as well as Federrath et al. 2010b). Our code also includes a sink
particle formation algorithm (Bate et al., 1995; Jappsen et al., 2005; Greif
et al., 2011) to properly track the formation of stars during each run.
Although our numerical models show active star formation by the end of the
simulations, we ignore the sink particles in this study, since we only want to
focus on the physical properties of the different gas phases. A more detailed
analysis of the star formation history of our runs and the sink particle
algorithm that we use can be found in Bertram et al. (2015).

We assume that the gas has a turbulent velocity field with an initial power
spectrum $P(k)\propto k^{-4}$ and consists of a natural mixture of solenoidal
and compressive modes (see, e.g. Federrath et al. 2010a). The strength of the
turbulence is set by our choice of the initial virial $\alpha$ parameter. This
is defined as $\alpha=E_{\text{kin}}/|E_{\text{pot}}|$, i.e. the kinetic
energy $E_{\text{kin}}=1/2M_{\text{tot}}\sigma_{v}^{2}$ divided by the
potential energy $E_{\text{pot}}=-3GM_{\text{tot}}^{2}/(5R)$ measured at the
beginning of each run. The quantity $\sigma_{v}$ denotes the 3D rms velocity
dispersion in the cloud. We note that the virial parameter is also often
calculated as $\alpha=2E_{\text{kin}}/|E_{\text{pot}}|$ in the literature,
which is different from the notation used in this paper (see, e.g. Federrath &
Klessen, 2012). With our definition, a value of $\alpha=0.5$ corresponds to a
cloud in virial equilibrium and $\alpha>1.0$ describes clouds that are
gravitationally unbound. In order to span a large range in the virial
parameter space, we analyze models with $\alpha=0.5,2.0\,$ and $8.0$. Note
that although it is convenient to parameterize the models in terms of
$\alpha$, the underlying physical quantity that changes as we change $\alpha$
is the turbulent velocity dispersion of the gas, since we keep $M_{\rm tot}$
and the initial radius $R$ fixed in all three models. It is therefore the
changes to the velocity field, and the consequent changes to the density
distribution of the gas, that dictate the observed trends seen in our
simulations. However, for simplicity and clarity, we will simply refer to the
different $\alpha$ models in our discussions in the following sections.

The turbulence is not driven and hence decays throughout the simulation (Mac
Low et al., 1998). The timesteps between the individual snapshots are $\Delta
t_{\text{snap}}\approx 7\,$kyr, corresponding to $\sim 200$ snapshots per
simulation in total. Table 1 summarizes the initial conditions of our clouds.
We note that due to the high computational cost of our $\alpha=0.5$ model, we
had to stop our run already at $t\approx 0.7\,t_{\text{ff}}$, while the other
runs end at one free-fall time. This is due to an intense refinement of dense
gas regions within the cloud. Furthermore, regarding the densities in the CMZ,
we have to keep in mind that a significant number of CMZ clouds have densities
higher than those modelled in our simulations, which would lead to even larger
turbulent velocities in our simulations for the same constant total mass.

###  2.3 Radiative transfer

We post-process our data in order to generate synthetic line emission maps for
the atomic and molecular tracers listed in Table 2, using the radiative
transfer code Radmc-3d111www.ita.uni-
heidelberg.de/$\sim$dullemond/software/radmc-3d/ (Dullemond, 2012). We use the
Large Velocity Gradient (LVG) approximation (Sobolev, 1957) to compute the C+,
O and CO level populations. The LVG implementation in Radmc-3d is described in
Shetty et al. (2011a). We use the rate coefficients for collisional excitation
and de-excitation by atomic and molecular hydrogen tabulated in the Leiden
database (Schöier et al., 2005).

At the moment, Radmc-3d cannot deal with Arepo data directly. Thus, we have to
map the simulation output onto a cubic grid. For that, we adopt a resolution
of 512 cells in all three spatial dimensions, corresponding to a cell size of
$\Delta x=\Delta y=\Delta z=0.087$ pc. Moreover, we also use 512 channels in
velocity space for our radiative transfer post-processing, corresponding to a
spectral resolution of 0.057 km s-1, 0.066 km s-1 and 0.102 km s-1 for our
models with $\alpha=0.5,2.0$ and 8.0. Our choice of spatial resolution is
determined by technical limitations within Radmc-3d: as it is a serial code,
our grid size is limited by the requirement that the data should all fit
within the memory of the computer that we use to carry our the post-
processing. However, for most of our tracers, we find only minor differences
in the emission as we increase the grid resolution. In particular, we find
that a resolution of $512^{3}$ grid cells is enough in order to properly
recover the emission of the cloud for all of the tracers apart from the [Oi]
63$\mu$m line (see Appendix A). Hence, our main conclusions are not strongly
affected by our choice of spatial resolution.

Our simulations do not explicitly track the abundance of 13CO and so we need a
procedure to relate the 13CO number density to that of 12CO. A common
assumption is that the ratio of 12CO to 13CO is identical to the elemental
abundance ratio of 12C to 13C (see, e.g. Roman-Duval et al., 2010). We make
the same assumption and set the 12CO to 13CO ratio to a constant value,
$R_{12/13}=60$. In reality, physical effects like chemical fractionation
(Watson et al., 1976) and photodissociation (Visser et al., 2009) will alter
the abundances of 13CO relative to 12CO, leading to a variable value of
$R_{12/13}$ within the cloud. However, Szűcs et al. (2014) have shown that the
resulting 13CO maps do not differ greatly from those that we would obtain by
assuming a constant ratio of 12CO to 13CO, and so we expect this approximation
to be adequate.

The radiative transfer calculation yields position-position-velocity (PPV)
cubes of brightness temperatures $T_{B}$, which are related to the intensity
via the Rayleigh-Jeans approximation,

$T_{B}(\nu)=\left(\frac{c}{\nu}\right)^{2}\frac{I_{\nu}}{2k_{B}},$ (2)

where $I_{\nu}$ is the specific intensity at frequency $\nu$ and $k_{B}$ the
Boltzmann constant. For each run, we carry out the post-processing for each of
the tracers listed in Table 2 along an arbitrary line of sight (LoS), using
the last snapshot of each model. Since we use isotropic turbulence without
magnetic fields, we do not expect our results to significantly depend on the
specific choice of the LoS and hence we only focus on emission observed along
the $z$-direction (see also Appendix B).

Tracer | Type | Transition | $\lambda$ | $\nu$
---|---|---|---|---
|  |  | [$\mu$m] | [GHz]
12CO | Molecular | $J=1\rightarrow 0$ | 2600 | 115
13CO | Molecular | $J=1\rightarrow 0$ | 2720 | 110
[Cii] | Atomic |  2$P$${}_{3/2}\rightarrow^{2}$$P$1/2 | 158 | 1900
[Oi] | Atomic |  3$P$${}_{1}\rightarrow^{3}$$P$2 | 63 | 4744
[Oi] | Atomic |  3$P$${}_{0}\rightarrow^{3}$$P$1 | 145 | 2060
Table 2: Tracer name and type, quantum mechanical transition, wavelength
$\lambda$ and frequency $\nu$ of various fine structure lines, which we model
with the radiative transfer code Radmc-3d.

##  3 Results

In this section, we analyze various physical cloud parameters using the latest
time snapshot (see Table 1). We start with the thermal state of the MCs
(Section 3.1) and investigate which chemical components best trace the total
shape of the cloud (Section 3.2). Afterwards, we estimate the effective cloud
radii by using synthetic observations in various diagnostic lines (Section
3.3). We then explore the ability of the lines listed in Table 2 to trace
dense MC regions (Section 3.4) and continue with a quantitative estimate of
the mass fraction of the cloud that is traced by the line emission (Section
3.5). At the end, we discuss the $X_{\text{CO}}$-factor of the MCs.

###  3.1 Thermal state of the clouds

In Fig. 1 we show the two-dimensional PDFs of gas temperatures and densities
within the clouds for the different virial parameters and for all chemical
components: total mass, molecular hydrogen, atomic oxygen, ionized carbon and
carbon monoxide. These PDFs are mass-weighted, with the colour-coding
indicating the fraction of the total mass of the cloud (in the upper row of
panels) or of the species of interest (in the remaining panels) located at
each point in density-temperature space. The number densities shown in the
upper row of panels are the total number density of particles; in the
remaining panels, we show instead the number density of the species of
interest ($n_{\rm H_{2}}$, $n_{\rm O}$, etc.). Figure 1 demonstrates that the
thermal state of the clouds for all of the components only weakly depends on
the virial parameter. However, we find strong differences in the PDFs between
the various tracers.

We find a significant amount of warm gas at densities $n\lesssim 1000\,$cm-3
and temperatures between $T\approx 300-10^{4}\,$K. We also note that there is
a substantial scatter in the temperature at every value of the density in this
regime. This gas corresponds to the warm and tenuous envelope in the outer
regions of the clouds with low extinction. It is heated by the strong external
ISRF via the photoelectric effect. For densities higher than $n\gtrsim
1000\,$cm-3, we find an inverse relationship between the temperature and the
density. The temperature decreases with increasing density, owing to the
efficient self-shielding of the dense gas from the external radiation field.
Such high densities only occur in the inner regions of the clouds with
temperatures of $T<100\,$K, which are cold enough to result in active star
formation (for further information about the star formation history of our
numerical simulations, we refer the reader to Bertram et al. 2015). These are
also the parts of the cloud in which we find molecular H2 and CO gas in our
simulations.

If we compare the thermal state of the total mass to the thermal state of the
other tracers, we find that only the atomic oxygen and ionized carbon show a
similar behavior over the whole range of temperatures. This is because these
components exist in the diffuse cloud regions as well as in the denser, more
shielded gas. We also see that molecular hydrogen is strongly dissociated in
the diffuse, warm envelope of the cloud. It only reaches maximal temperatures
of $\sim 4000\,$K, which is lower than the warmest regions observed for the
total gas. In the case of CO, the range of temperatures ($30-300\,$K) is even
smaller. This suggests that CO traces a denser, colder part of the cloud than
[Cii] or [Oi]. We will examine the consequences of this in more detail in the
following sections.

Figure 1: Physical state of the clouds at the end of our simulations. We show
two-dimensional PDFs of gas temperature and number density for various
chemical components for runs with $\alpha=0.5$ (left column), 2.0 (middle
column), and 8.0 (right column). In the upper row, the number densities shown
are the total particle number density. In the remaining rows, on the other
hand, they are the number density of the chemical species indicated in the
panel (e.g. $n_{\rm H_{2}}$, $n_{\rm O}$, etc.). Note that the scaling of the
density axis therefore differs from panel to panel. The color coding shows the
fraction of the total gas mass in each logarithmic density and temperature
bin.

###  3.2 Tracing the shape of the MC

Figure 2: Logarithmic column density maps computed along the LoS in the
$z$-direction for the total density (top row) and the H2 density (bottom row)
for different virial parameters: $\alpha=0.5$ (left column), 2.0 (middle
column) and 8.0 (right column). Each side has a length of $44.5\,$pc. The
initial cloud radius of 8.9 pc is indicated in the bottom panels.

Figure 3: Total mass per unit logarithmic column density for the total and the H2 gas as shown in Fig. 2 for our different virial parameters: $\alpha=0.5$ (top), 2.0 (middle) and 8.0 (bottom). Model | [Cii] (158 $\mu$m) | [Oi] (145 $\mu$m) | [Oi] (63 $\mu$m) |  12CO (2600 $\mu$m) |  13CO (2720 $\mu$m)
---|---|---|---|---|---
| [K km s-1] | [K km s-1] | [K km s-1] | [K km s-1] | [K km s-1]
GC-0.5-1000 | 3.0 | 1.0 | 0.6 | 1.0 | 1.0
GC-2.0-1000 | 3.0 | 1.0 | 0.6 | 0.3 | 0.3
GC-8.0-1000 | 3.0 | 0.8 | 0.5 | 0.1 | 0.1
Table 3: Threshold values $W_{\text{thresh}}$ used for the integrated intensity maps in Fig. 4 in order to estimate the effective cloud radius $R_{\text{eff}}$. The values are chosen such that the denser cloud regions can be clearly separated from the diffuse gas phase around it, leading to a reasonable estimate for the effective cloud radii presented in Table 4. Model | 80% of total mass | [Cii] (158 $\mu$m) | [Oi] (145 $\mu$m) | [Oi] (63 $\mu$m) |  12CO (2600 $\mu$m) |  13CO (2720 $\mu$m)
---|---|---|---|---|---|---
| [pc] | [pc] | [pc] | [pc] | [pc] | [pc]
GC-0.5-1000 | 11.0 | 9.0 | 8.7 | 8.1 | 6.6 | 6.2
GC-2.0-1000 | 13.1 | 13.8 | 13.0 | 12.2 | 7.6 | 6.8
GC-8.0-1000 | 17.5 | 20.1 | 19.9 | 19.7 | 8.4 | 6.9
Table 4: Effective radii of the clouds shown in Fig. 4, as inferred from
equation 3, based on the threshold intensities given in Table 3. For
comparison, we also give the effective radii of the clouds that incorporate
$\sim 80$% of the total gas mass in the second column.

Fig. 2 shows logarithmic column density maps for the total density and the H2
density computed along the LoS in the $z$-direction for virial parameters
$\alpha=0.5,2.0$ and 8.0 for the last time snapshots (see Table 1). We see
that as we increase $\alpha$, and hence the velocity dispersion of the gas,
the structure of the cloud changes significantly. When $\alpha=0.5$ and the
cloud is gravitationally bound, it remains relatively compact. It is therefore
able to shield itself relatively well from the effects of the high external
ISRF, with the result that the H2 column density traces the total column
density fairly well. In the runs with $\alpha=2$ and $\alpha=8$, however, the
cloud is not gravitationally bound and hence starts to expand. In addition,
the strong turbulence opens up large channels into the cloud by lowering the
column density locally, allowing radiation to penetrate deeper into the MC,
leading to enhanced photodissociation of the molecular hydrogen. Hence, we do
not find a good correlation between the total and the H2 column density in the
$\alpha=2.0$ and $\alpha=8.0$ runs. Molecular hydrogen can only exist in the
most dense regions of the cloud, where it is best able to self-shield from the
strong external radiation field.

This is also illustrated in Fig. 3, which shows the corresponding column
density PDFs for three virial parameter models. If we compare the H2 column
density PDFs with each other, we find that the peak column density in the
$\alpha=0.5$ model is shifted to larger values compared to the $\alpha=8.0$
model by about one order of magnitude. This is because the turbulent kinetic
energy in the virialized $\alpha=0.5$ run is much smaller than in the
$\alpha=8.0$ run, leading to denser cloud regions and thus to larger average
column densities. Moreover, the fraction of H2 gas in the total box decreases
with increasing virial parameter, owing to the high external ISRF, which
photodissociates those molecular components more effectively as the radiation
can penetrate deeper into the cloud’s interior. Furthermore, we see that the
H2 gas generally extends to much lower column densities than the total gas.
This is because the total gas consists of both atomic and molecular hydrogen
and thus does not fall below lower values than $N_{\text{tot}}\approx
10^{21}\,$cm-2 due to the Hi background, while the fraction of H2 gas can
continuously decrease in the outer regions due to photodissociation.

In analogy to Fig. 2, Fig. 4 shows velocity-integrated intensity maps computed
along the LoS in the $z$-direction for the same virial parameters and the
different tracers presented in Table 2. We again find that the molecular
tracers 12CO and 13CO are photodissociated at the edges of the MC by the
strong external ISRF. This effect is even stronger if the amount of turbulent
kinetic energy in the box is high. Comparing the different tracers to the
total gas column in Fig. 2, we find that the cloud’s shape is well reproduced
by all atomic tracers, although the integrated brightness temperatures of
these lines are significantly smaller than those of the 12CO and 13CO lines.
This holds for all three virial parameter models. Conversely, in the case of
CO, we find that neither of its isotopologues describes the total gas
distribution of the cloud. However, we observe that they do remain reasonably
good tracers of the dense H2, owing to the effect of self-shielding. We will
quantify this finding in Section 3.5.

###  3.3 Estimating the effective cloud radii

In the previous section we established that carbon monoxide is not a good
tracer for the total column density, while the atomic species better reflect
the distribution of the total gas mass. In this context, CO tends to
significantly underestimate the radius of the total cloud in such an extreme
environment. This can compromise any observational estimates of the virial
parameter $\alpha$, which requires an accurate estimate for the cloud’s radius
(see also our follow-up study in Paper II). However, we note that it is
generally complicated to define a proper cloud radius. Molecular clouds are
complex hierarchical systems and estimates of cloud radii always come with
some assumptions of how to define outer cloud boundaries. Hence, if we refer
to the effective cloud radius in the following, we always mean the radius of a
spherical cloud with the same surface area.

Below, we try to estimate the effective cloud radius $R_{\text{eff}}$ in a
similar fashion to what would be applied to observational data sets. At first
sight, this seems to be complicated, since all clouds shown in Fig. 4 are far
away from a spherical geometry. Nevertheless, we can estimate an effective
radius by adopting an integrated intensity threshold $W_{\text{thresh}}$ for
each map and counting the number of cells $N$, which fulfill
$W_{\text{cloud}}>W_{\text{thresh}}$. The values for $W_{\text{thresh}}$ are
listed in Table 3 and are chosen so that there is a clear separation between
emission from the cloud and from the diffuse gas in which it is embedded. Our
results are not particularly sensitive to the exact value chosen for
$W_{\text{thresh}}$, as long as it is large enough to not be contaminated by
the diffuse emission, but not so large that we miss significant emission from
the cloud itself. For example, in the case of model GC-0.5-1000, values of
$W_{\text{thresh}}$ anywhere in the range $2<W_{\text{thresh}}<10\,$K km s-1
give essentially the same results for $R_{\text{eff}}$. For comparison, we
also give the effective radii of the clouds that incorporate $\sim 80$% of the
total gas mass.

We have also verified that our adapted thresholds are comparable or greater
than the practical thresholds that we would obtain using modern instruments to
measure the line emission. For example, in a single pointing, the GREAT
instrument (German REceiver for Astronomy at Terahertz Frequencies, Heyminck
et al., 2012) on-board the Stratospheric Observatory for Infrared Astronomy
(SOFIA) can measure integrated intensities of the [Cii] line down to around
$\sim 0.06\,$K km s-1 in around one hour of integration
time222https://great.sofia.usra.edu/cgi-bin/great/great.cgi, although
multipixel array detectors are available today that become more and more
efficient. Mapping an entire cloud using multiple pointings obviously allows
one to spend much less time per pointing, yielding a higher threshold, but it
remains plausible to map extended regions with a sensitivity comparable to or
better than the adopted threshold.

We compute the effective cloud radius $R_{\text{eff}}$ by assuming an
effective cloud area $A_{\text{eff}}=\pi R_{\text{eff}}^{2}$, where
$A_{\text{eff}}=N\Delta x^{2}$ and $\Delta x$ is the cell size introduced in
Section 2.3. Thus, we get

$R_{\text{eff}}=\sqrt{\frac{A_{\text{eff}}}{\pi}}=\Delta
x\sqrt{\frac{N}{\pi}}.$ (3)

Table 4 shows the different effective radii of the various clouds presented in
Fig. 4, computed via equation (3). If we compare the different values listed
in Table 4, we find that the cloud radii measured with the integrated
intensities of 12CO and 13CO are significantly smaller by a factor of
$1.5-2.5$ compared to the values measured with the integrated intensities of
the atomic tracers. This is because the strong ISRF and high CRF destroy the
molecular tracers in the more diffuse regions, leading to smaller effective
radii. Conversely, the radii measured with the atomic components roughly agree
with each other for one specific virial parameter model and give the best
radius estimate for the total cloud that incorporates $\sim 80$% of all gas.
These values confirm our assumption that the atomic species accurately reflect
the spatial distribution of the total gas column and that they yield a robust
estimate of the effective radius of the cloud.

###  3.4 Tracing the total and H2 column density

Figure 4: Velocity-integrated intensity maps computed along the LoS in the
$z$-direction for different virial parameters: $\alpha=0.5$ (left column), 2.0
(middle column) and 8.0 (right column). From top to bottom: integrated
intensity maps for 12CO ($J=1\rightarrow 0$), 13CO ($J=1\rightarrow 0$), [Oi]
(63 $\mu$m), [Oi] (145 $\mu$m) and [Cii] (158 $\mu$m). Note the different
scaling in the colorbars at the right hand side. Each side has a length of
$44.5\,$pc. The initial cloud radius of 8.9 pc is indicated in the bottom
panels. The contour lines show the threshold values given in Table 3.

Figure 5: Log-log plot showing velocity-integrated intensities against total
column density for different virial parameters: $\alpha=0.5$ (left column),
2.0 (middle column) and 8.0 (right column). Shown are all our tracers
presented in Table 2. In each plot, the mean value is indicated by a black
solid line.

Figure 6: Same as Fig. 5, but with the column density of H2.

Now we turn to the densest parts of the cloud, because this is where stars are
formed, and we study how well atomic fine-structure lines can be used to trace
this regime. Fig. 5 and 6 show velocity-integrated intensities plotted against
the total column density and H2 column density for the different virial
parameters and all tracers presented in Table 2. In these plots, we also show
the the mean value at each column density, which is indicated by a black solid
line. The small fluctuations at high column densities in each plot are caused
by the small number of cells available to compute an average value of $W$.

We find clear differences for the various chemical components. For the total
density in Fig. 5, we observe that the 12CO and 13CO tracers fall off sharply
with decreasing column density below $\sim 10^{22}\,$cm-2, which is due to
photodissociation of the gas in the diffuse cloud regions, owing to the strong
external radiation field. On the other hand, at high column densities, we find
a saturation of the 12CO emission due to the effect of the line opacity, while
the 13CO emission is a significantly better tracer of the column density than
12CO. This is because 13CO is optically thin and consequently can better trace
compact regions within the cloud. However, even in a small regime where
$N_{\text{tot}}$ and $W_{\text{CO}}$ show some degree of correlation, we find
a substantial amount of scatter around the mean value for both 12CO and 13CO.
Furthermore, we also show the scatter plots for the different atomic species
in Fig. 5. As expected, we see that all atomic components extend to
significantly lower column densities of $\sim 10^{21}\,$cm-2 compared to
carbon monoxide. Hence, both atomic carbon and oxygen are much better tracers
of the low column density material of the total gas in the cloud. However,
both the [Cii] 158 $\mu$m and the [Oi] 63 $\mu$m line saturate at higher
column densities, owing to the effect of the line opacity and thus either 13CO
or the [Oi] 145 $\mu$m line is to be preferred as tracers for the dense cloud
material.

For the H2 column density in Fig. 6, we find a different result. Firstly, we
see that both 12CO and 13CO extend to significantly lower H2 column densities
of $\sim 5\times 10^{20}-10^{21}\,$cm-2, reflecting the fact that CO is a
better tracer of molecular mass than of total mass. At higher column
densities, we find a qualitatively similar behavior as for the total gas in
Fig. 5, again owing to the effect of the line opacities. Moreover, we also
show the scatter plots for the different atomic species. We note that both the
atomic oxygen and carbon extend over a much larger range in H2 column density
than in total gas column density, showing that those components can also be
used in order to accurately trace a large fraction of the molecular hydrogen
mass. However, the correlation between the emission of [Cii] 158 $\mu$m, [Oi]
63 $\mu$m and $N_{\text{H}_{2}}$ breaks down at H2 column densities above
$\sim 10^{20}-10^{21}\,$cm-2, owing to the optical depth of these lines in
this density regime. This means that one has to rely on other probes in order
to infer information about the compact gas regions. On the contrary, the
correlation between the [Oi] 145 $\mu$m emission and $N_{\text{H}_{2}}$ seems
to be provide a better way to study the dense H2 clumps in the cloud.

We also analyze the impact that our over-estimation of the CO abundance may
have on the relationship between the integrated intensity of the atomic
tracers and the column density (see Appendix C). However, the influence of
this turns out to be small.

###  3.5 Mass fraction traced by the emission

Figure 7: Cumulative fraction of the total (top row) and H2 (bottom row) mass
in the cloud for different virial parameters $\alpha=0.5$ (left column), 2.0
(middle column) and 8.0 (right column), traced by our various chemical
components presented in Table 2, as a function of the minimum velocity-
integrated intensity in the line.

In Section 3.2 we speculated that a significant fraction of the total gas mass
may not be traced by 12CO or 13CO. To quantify this finding, Fig. 7 shows the
fraction of total and H2 gas mass observed along the LoS having velocity-
integrated intensities greater than a minimum value $W_{\text{min}}$ for our
various tracers. Furthermore, we analyze how this mass fraction changes as we
vary $W_{\text{min}}$. In this context, Fig. 7 defines a theoretical limit for
detecting gas above a given sensitivity threshold, e.g. determined by the
telescope or the detector.

We see that all lines are a strong function of the minimum velocity-integrated
intensity and of the virial parameter. In case of the mass fraction of the
total gas, we find that our different atomic tracers recover almost all of the
total mass if the minimum integrated intensity is $W_{\text{min}}\approx
0.1\,$K km s-1 in the $\alpha=0.5$ model and $W_{\text{min}}\approx 1.0\,$K km
s-1 in the $\alpha=2.0$ and $\alpha=8.0$ models. Furthermore, as already
indicated in Section 3.2, we see that in the case of very small threshold
values $W_{\text{min}}$ in the $\alpha=0.5$ model, both 12CO and 13CO only
trace about 75% of the total mass or less. This fraction decreases even
further for higher $\alpha$ values. For example, if we adopt
$W_{\text{min}}\approx 0.1\,$K km s-1 in the cloud with $\alpha=8.0$, we find
that only $\sim 40\%$ of the total mass is traced by CO, while all atomic
tracers recover almost 100% of the total gas mass. Furthermore, we observe
narrow emission ranges for some of the atomic components, which trace almost
100% of the total mass. For example, regarding our intermediate $\alpha=2.0$
model and the [Oi] 63 $\mu$m tracer, almost 100% of the total mass lies in the
narrow emission range between $0.5<W_{\text{min}}<3.0$. In the case of the
other $\alpha$ models, we find a similar behavior. This is because there are
almost no lines-of-sight having higher values of integrated intensity (see
also Fig. 5 and 6), which may be different in even denser clouds.

Regarding the H2 mass, we see that both 12CO and 13CO trace a higher H2 mass
fraction compared to the total mass of the cloud. This is because the
molecular hydrogen can only exist in dense regions of the cloud, owing to the
ability to self-shield from the external radiation field. This is similar to
our finding for the carbon monoxide. Furthermore, we observe that the atomic
components recover almost 100% of the H2 mass if we set $W_{\text{min}}\approx
0.5\,$K km s-1. Hence, we find that those components again trace a
significantly higher mass fraction of the gas compared to the carbon monoxide.
Regarding both 12CO and 13CO, we see that the amount of mass recovered by
those molecular tracers decreases with higher virial parameter for a fixed
minimum value of the integrated intensity. For example, if we adopt a value of
$W_{\text{min}}\approx 0.1\,$K km s-1, almost 98% of the H2 mass is traced by
12CO for our $\alpha=0.5$ model, while $\sim 90\%$ is traced in our
$\alpha=2.0$ and only $\sim 70\%$ can be recovered in the $\alpha=8.0$ model.
Moreover, similar to the mass fraction of the total gas mass, we also find a
narrow emission range for the various atomic lines in all $\alpha$ models,
which trace about 100% of the H2 mass. In contrast, the emission of the 12CO
and 13CO molecules is extended over a much wider range in integrated
intensities. Hence, a much larger emission range is needed in order to trace
the whole mass of molecular hydrogen in the cloud.

In summary, we conclude that CO is only a good tracer for the molecular
content of the cloud, but can miss a significant fraction of the total gas
mass, particularly for unbound clouds.

###  3.6 Estimating the CO-to-H2 conversion factor

Figure 8: Logarithmic maps of the CO-to-H2 conversion factor,
$X_{\text{CO}}=N_{\text{H}_{2}}/W_{\text{CO}}$, for the three models:
$\alpha=0.5$ (top), 2.0 (middle) and 8.0 (bottom). The average
$X_{\text{CO}}$-factor values are given in Table 5. We see that only the inner
parts of the cloud yield the canonical $X_{\text{CO}}$-factor value on
average, while the outer cloud regions exhibit significantly larger
$X_{\text{CO}}$ values. The gray background denotes regions where the CO
emission is zero. Hence, it is impossible to compute a value for
$X_{\text{CO}}$ there.

Based on our maps of column densities for H2 (Fig. 2) and velocity-integrated
intensities for 12CO (Fig. 4), we also estimate the value of the CO-to-H2
conversion factor, $X_{\text{CO}}$, in our various models. $X_{\text{CO}}$ is
a widely used quantity in particular in extragalactic astronomy to derive H2
column densities from CO observations (Bolatto et al., 2013). It is defined
via

$X_{\text{CO}}=\frac{N_{\text{H}_{2}}}{W_{\text{CO,1-0}}},$ (4)

where $N_{\text{H}_{2}}$ denotes the H2 column density and $W_{\text{CO,1-0}}$
the velocity-integrated brightness-temperature. Previous studies tried to
estimate the $X_{\text{CO}}$-factor in both observations and numerical
simulations (see, e.g. Solomon et al., 1987; Polk et al., 1988; Young &
Scoville, 1991; Dame et al., 2001; Liszt et al., 2010; Leroy et al., 2011;
Glover & Mac Low, 2011; Glover & Clark, 2012; Bolatto et al., 2013; Narayanan
& Hopkins, 2013; Clark & Glover, 2015). For example, Shetty et al. (2011a,b)
investigated the $X_{\text{CO}}$-factor in numerical simulations of turbulent
clouds, varying environmental properties such as the initial number density,
the metallicity and the external UV field, finding cloud average values
$X_{\text{CO}}\approx 2-4\times 10^{20}\,$cm-2 K-1 km-1 s for solar-
metallicity models. Moreover, Bolatto et al. (2013) review the efforts to
measure the CO-to-H2 conversion factor for molecular clouds in the Milky Way
disk, recommending a similar factor $X_{\text{CO}}\sim 2\times 10^{20}\,$cm-2
K-1 km-1 s with $\pm 30\%$ uncertainty. However, all these studies arrived at
the conclusion that the $X_{\text{CO}}$-factor in the Galaxy is remarkably
constant, having a value of $X_{\text{CO}}\approx 2-4\times 10^{20}\,$cm-2 K-1
km-1 s.

Nevertheless, in the last few years concerns have been raised that the
$X_{\text{CO}}$-factor may be significantly different from the canonical value
in MCs located in the galactic center. For example, Narayanan et al. (2012)
have shown that large gas temperatures and velocity dispersions can increase
the CO intensity, thus decreasing $X_{\text{CO}}$. On the other hand, the
strong ISRF and CRF can also destroy the carbon monoxide, increasing
$X_{\text{CO}}$ (Clark & Glover, 2015). In addition, observational
measurements of $X_{\text{CO}}$ at the center of other nearby spiral galaxies
typically find values that are lower than the canonical Galactic value
(Sandstrom et al., 2013). Moreover, conditions in the Galactic Center are
somewhat similar to those found in ULIRGs, and there is considerable
observational evidence that $X_{\text{CO}}$ is smaller in ULIRGs than in
normal spiral galaxies (see e.g. the detailed discussion in Bolatto et al.,
2013). So far, we are still missing a reliable picture of the
$X_{\text{CO}}$-factor in different (and extreme) physical environments.

Model | $X_{\text{CO}}$
---|---
| [cm-2 K-1 km-1 s]
GC-0.5-1000 | $3.9\times 10^{20}$
GC-2.0-1000 | $2.6\times 10^{20}$
GC-8.0-1000 | $1.3\times 10^{20}$
Table 5: Values of the CO-to-H2 conversion factor, relating the amount of 12CO
($J=1\rightarrow 0$) emission $W_{\text{CO}}$ to the H2 column density
$N_{\text{H}_{2}}$ for all our $\alpha$ runs. For the computation of the
$X_{\text{CO}}$-factor, we only include pixels where the velocity-integrated
intensity for carbon monoxide is larger than the corresponding threshold given
in Table 3. Thus, $X_{\text{CO}}$ only accounts for regions with a significant
amount of CO emission in the cloud.

In the following, we evaluate the $X_{\text{CO}}$-factor for our three
$\alpha$ models and list our results in Table 5. In general, the
$X_{\text{CO}}$-factors given in Table 5 are computed by taking the ratio of
the mean values of the H2 column density and the velocity-integrated intensity
of CO. We only include pixels where the velocity-integrated intensity for
carbon monoxide is larger than the corresponding threshold given in Table 3.
These are 1.0, 0.3 and 0.1 K km s-1 for the models with $\alpha=0.5,2.0$ and
$8.0$. Thus, $X_{\text{CO}}$ only accounts for regions with a significant
available amount of CO emission in the cloud, as illustrated in Fig. 4.

We find a range of values for the CO-to-H2 conversion factor in Table 5,
reaching from $\sim 1-4\times 10^{20}\,$cm-2 K-1 km-1 s. Our estimates are in
good agreement to estimates of the canonical $X_{\text{CO}}$-factor in the
Milky Way. However, we again emphasize that our Nelson & Langer (1997) network
overestimates the rate at which CO forms (see also Section 2.1) and hence we
expect that our $X_{\text{CO}}$-factors listed in Table 5 are lower limits.

Furthermore, we show logarithmic maps of the $X_{\text{CO}}$-factor for all
our $\alpha$ models in Fig. 8. Here, the $X_{\text{CO}}$-values are evaluated
by computing the ratio of the H2 column density and the velocity-integrated
intensity of CO for each individual pixel. The gray background denotes regions
where the CO emission is zero. Hence, it is impossible to compute a value for
$X_{\text{CO}}$ there. We find significantly different values for the inner
and the outer parts of the clouds. While the inner cloud regions reproduce the
canonical $X_{\text{CO}}$-factor value on average, we find that the outer
cloud regions exhibit significantly larger $X_{\text{CO}}$ values compared to
the canonical one. This is because the carbon monoxide is photodissociated at
the edges of our clouds by the strong ISRF, leading to an increase in
$X_{\text{CO}}$. In contrast, the molecular CO is better able to self-shield
in the dense parts of the cloud and hence we observe that those regions fairly
reproduce the canonical value.

##  4 Summary and Conclusions

In this paper we have analyzed synthetic images of MCs in a CMZ-like
environment. For this purpose, we have performed numerical simulations of
model clouds with the moving mesh code Arepo (Springel, 2010) using
environmental properties comparable to those experienced by typical CMZ
clouds. We adopted values for the interstellar radiation field (ISRF) and the
cosmic-ray flux (CRF) that are a factor of $\sim 1000$ larger than the values
measured in the solar neighbourhood (Clark et al., 2013). We simulated clouds
with an initial number density of $n_{0}=10^{3}\,$cm-3 and studied the impact
of different virial $\alpha$ parameters of $\alpha=0.5,2.0$ and $8.0$. The
total cloud mass was set to a constant value of $M_{\text{tot}}=1.3\times
10^{5}\,$M⊙. Furthermore, we used the radiative transfer code Radmc-3d
(Dullemond, 2012) to compute synthetic maps of our clouds in important
diagnostic lines, i.e. we model the emission of the cloud in [Cii] (158
$\mu$m), [Oi] (145 $\mu$m), [Oi] (63 $\mu$m), 12CO (2600 $\mu$m) and 13CO
(2720 $\mu$m). We report the following findings:

  * •

Atomic carbon and oxygen is found in both cold and warm regions in the clouds,
but CO only traces cold regions (see Section 3.1).

  * •

If the cloud is virialized ($\alpha=0.5$), the H2 gas is much better able to
self-shield from the high external ISRF compared to other runs with a larger
$\alpha$ value. This is because the gravitational pressure compresses the gas
to much higher densities because the amount of turbulent kinetic energy in the
box is smaller in the low $\alpha$ models than in the high virial parameter
models. In the latter case, radiation can penetrate deeper into the MC,
leading to photodissociation of the molecular hydrogen (see Section 3.2).

  * •

We show that the atomic components trace the shape of the total cloud very
well even in diffuse cloud regions (see Section 3.2).

  * •

On the other hand, the molecular gas (H2 as well as CO) is photodissociated at
the edges of the MC by the strong external ISRF.

  * •

As a consequence, the photodissociation has a strong impact on measurements of
the effective MC radius using a molecular tracer, which tend to significantly
underestimate the radius of the total cloud in environments where the external
radiation field and the turbulent kinetic energy is high. However, we find
that the atomic tracers recover a significantly larger fraction of the total
cloud mass (see Section 3.3).

  * •

We find a large dynamical range for the [Oi] 145 $\mu$m line in column density
space, which traces both the H2 and the total gas up to the most dense regions
in the cloud. In contrast, the [Oi] 63 $\mu$m and [Cii] 158 $\mu$m lines
already saturate at significantly lower column densities. Hence, we would not
expect that those components trace the compact MC regions very well (see
Section 3.4).

  * •

We find that all atomic components trace almost 100% of the H2 and the total
gas mass above a sensitivity threshold of $W_{\text{min}}\approx 0.1\,$K km
s-1 for the velocity-integrated intensities. However, we also find that CO
only traces a significantly lower mass fraction of H2 and the total density
compared to the atomic components (see Section 3.5).

  * •

We compute values of the CO-to-H2 conversion factor
$X_{\text{CO}}=N_{\text{H}_{2}}/W_{\text{CO}}$, which relates the H2 column
density $N_{\text{H}_{2}}$ to the amount of 12CO ($J=1\rightarrow 0$) emission
$W_{\text{CO}}$. We find values in the range $X_{\text{CO}}\approx 1-4\times
10^{20}\,$cm-2 K-1 km-1 s, in good agreement to the canonical
$X_{\text{CO}}$-factor values obtained in observations of MCs in the Milky Way
(see Section 3.6). However, we caution that simplifications in our chemical
model mean that these values are lower limits.

## Acknowledgements

We thank the referee for an insightful and very constructive report. We also
thank Stella Offner and Diane Cormier for informative discussions about the
project. EB, SCOG and RSK acknowledge support from the Deutsche
Forschungsgemeinschaft (DFG) via the SFB 881 (sub-projects B1, B2 and B8) “The
Milky Way System”, and the SPP (priority program) 1573, “Physics of the ISM”.
Furthermore, EB acknowledges financial support from the Konrad-Adenauer-
Stiftung (KAS) via their “Promotionsförderung”. The simulations presented in
this paper were performed on the Milkyway supercomputer at the Jülich
Forschungszentrum, funded via SFB 881. The radiative transfer post-processing
was performed on the kolob cluster at the University of Heidelberg, which is
funded in part by the DFG via Emmy-Noether grant BA 3706, and via a Frontier
grant of Heidelberg University, sponsored by the German Excellence Initiative
as well as the Baden-Württemberg Foundation. RSK acknowledges support from the
European Research Council under the European Community’s Seventh Framework
Programme (FP7/2007-2013) via the ERC Advanced Grant ”STARLIGHT: Formation of
the First Stars” (project number 339177). SER acknowledges support from
VIALACTEA, a Collaborative Project under Framework Programme 7 of the European
Union, funded under Contract #607380.

## References

  * Bania (1977) Bania T. M., 1977, ApJ, 216, 381
  * Bate et al. (1995) Bate M. R., Bonnell I. A., Price N. M., 1995, MNRAS, 277, 362
  * Bertram et al. (2015) Bertram E., Glover S. C. O., Clark P. C., Klessen R. S., 2015, MNRAS, 451, 3679
  * Bisbas et al. (2015) Bisbas T. G., Papadopoulos P. P., Viti S., 2015, ApJ, 803, 37
  * Bitran et al. (1997) Bitran M., Alvarez H., Bronfman L., May J., Thaddeus P., 1997, A&AS, 125, 99
  * Bolatto et al. (2013) Bolatto A. D., Wolfire M., Leroy A. K., 2013, ARA&A, 51, 207
  * Burton et al. (1978) Burton W. B., Liszt H. S., Gordon M. A., 1978, J. Phys. F, 39, 1
  * Clark & Glover (2015) Clark P. C., Glover S. C. O., 2015, MNRAS, 452, 2057
  * Clark et al. (2012) Clark P. C., Glover S. C. O., Klessen R. S., 2012, MNRAS, 420, 745
  * Clark et al. (2013) Clark P. C., Glover S. C. O., Ragan S. E., Shetty R., Klessen R. S., 2013, ApJ, 768, L34
  * Dame et al. (2001) Dame T. M., Hartmann D., Thaddeus P., 2001, ApJ, 547, 792
  * Draine (1978) Draine B. T., 1978, ApJS, 36, 595
  * Dullemond (2012) Dullemond C. P., , 2012, RADMC-3D: A multi-purpose radiative transfer tool, Astrophysics Source Code Library
  * Federrath et al. (010b) Federrath C., Banerjee R., Clark P. C., Klessen R. S., 2010b, ApJ, 713, 269
  * Federrath & Klessen (2012) Federrath C., Klessen R. S., 2012, ApJ, 761, 156
  * Federrath et al. (010a) Federrath C., Roman-Duval J., Klessen R. S., Schmidt W., Mac Low M.-M., 2010a, A&A, 512, A81
  * Glover & Clark (2012) Glover S. C. O., Clark P. C., 2012, MNRAS, 421, 116
  * Glover et al. (2010) Glover S. C. O., Federrath C., Mac Low M.-M., Klessen R. S., 2010, MNRAS, 404, 2
  * Glover & Mac Low (2007) Glover S. C. O., Mac Low M.-M., 2007, ApJS, 169, 239
  * Glover & Mac Low (2011) Glover S. C. O., Mac Low M.-M., 2011, MNRAS, 412, 337
  * Greif et al. (2011) Greif T. H., Springel V., White S. D. M., Glover S. C. O., Clark P. C., Smith R. J., Klessen R. S., Bromm V., 2011, ApJ, 737, 75
  * Güsten et al. (1981) Güsten R., Walmsley C. M., Pauls T., 1981, A&A, 103, 197
  * Habing (1968) Habing H. J., 1968, Bulletin of the Astronomical Institutes of the Netherlands, 19, 421
  * Heyminck et al. (2012) Heyminck S., Graf U. U., Güsten R., Stutzki J., Hübers H. W., Hartogh P., 2012, A&A, 542, L1
  * Immer et al. (2012) Immer K., Menten K. M., Schuller F., Lis D. C., 2012, A&A, 548, A120
  * Indriolo et al. (2015) Indriolo N., Neufeld D. A., Gerin M., Schilke P., Benz A. O., Winkel B., Menten K. M., Chambers E. T., Black J. H., Bruderer S., et al. 2015, ApJ, 800, 40
  * Jappsen et al. (2005) Jappsen A.-K., Klessen R. S., Larson R. B., Li Y., Mac Low M.-M., 2005, A&A, 435, 611
  * Johnston et al. (2014) Johnston K. G., Beuther H., Linz H., Schmiedeke A., Ragan S. E., Henning T., 2014, A&A, 568, A56
  * Kauffmann et al. (2013) Kauffmann J., Pillai T., Zhang Q., 2013, ApJ, 765, L35
  * Leroy et al. (2011) Leroy A. K., Bolatto A., Gordon K., Sandstrom K., Gratier P., Rosolowsky E., Engelbracht C. W., Mizuno N., Corbelli E., Fukui Y., Kawamura A., 2011, ApJ, 737, 12
  * Lis & Menten (1998) Lis D. C., Menten K. M., 1998, ApJ, 507, 794
  * Lis et al. (1994) Lis D. C., Menten K. M., Serabyn E., Zylka R., 1994, ApJ, 423, L39
  * Lis et al. (2001) Lis D. C., Serabyn E., Zylka R., Li Y., 2001, ApJ, 550, 761
  * Liszt & Burton (1978) Liszt H. S., Burton W. B., 1978, ApJ, 226, 790
  * Liszt et al. (2010) Liszt H. S., Pety J., Lucas R., 2010, A&A, 518, A45
  * Longmore et al. (2013) Longmore S. N., Bally J., Testi L., Purcell C. R., Walsh A. J., Bressert E., Pestalozzi M., Molinari S., Ott J., Cortese L., Battersby C., Murray N., Lee E., Kruijssen J. M. D., Schisano E., Elia D., 2013, MNRAS, 429, 987
  * Longmore et al. (2012) Longmore S. N., Rathborne J., Bastian N., Alves J., Ascenso J., Bally J., Testi L., Longmore A., Battersby C., Bressert E., Purcell C., Walsh A., Jackson J., Foster J., Molinari S., Meingast S., Amorim A., Lima J., 2012, ApJ, 746, 117
  * Mac Low et al. (1998) Mac Low M.-M., Klessen R. S., Burkert A., Smith M. D., 1998, Physical Review Letters, 80, 2754
  * Martin et al. (2004) Martin C. L., Walsh W. M., Xiao K., Lane A. P., Walker C. K., Stark A. A., 2004, ApJS, 150, 239
  * Mathis et al. (1983) Mathis J. S., Mezger P. G., Panagia N., 1983, A&A, 128, 212
  * Molinari et al. (2011) Molinari S., Bally J., Noriega-Crespo A., Compiègne M., Bernard J. P., Paradis D., Martin P., Testi L., Barlow M., Moore T., Plume R., Swinyard B., Zavagno A., Calzoletti L., Di Giorgio A. M., Elia 2011, ApJ, 735, L33
  * Morris & Serabyn (1996) Morris M., Serabyn E., 1996, ARA&A, 34, 645
  * Narayanan & Hopkins (2013) Narayanan D., Hopkins P. F., 2013, MNRAS, 433, 1223
  * Narayanan et al. (2012) Narayanan D., Krumholz M. R., Ostriker E. C., Hernquist L., 2012, MNRAS, 421, 3127
  * Nelson & Langer (1997) Nelson R. P., Langer W. D., 1997, ApJ, 482, 796
  * Oka et al. (1998) Oka T., Hasegawa T., Sato F., Tsuboi M., Miyazaki A., 1998, ApJS, 118, 455
  * Polk et al. (1988) Polk K. S., Knapp G. R., Stark A. A., Wilson R. W., 1988, ApJ, 332, 432
  * Rathborne et al. (2015) Rathborne J. M., Longmore S. N., Jackson J. M., Alves J. F., Bally J., Bastian N., Contreras Y., Foster J. B., Garay G., Kruijssen J. M. D., Testi L., Walsh A. J., 2015, ApJ, 802, 125
  * Rathborne et al. (2014) Rathborne J. M., Longmore S. N., Jackson J. M., Foster J. B., Contreras Y., Garay G., Testi L., Alves J. F., Bally J., Bastian N., Kruijssen J. M. D., Bressert E., 2014, ApJ, 786, 140
  * Rodríguez-Fernández et al. (2004) Rodríguez-Fernández N. J., Martín-Pintado J., Fuente A., Wilson T. L., 2004, A&A, 427, 217
  * Roman-Duval et al. (2010) Roman-Duval J., Jackson J. M., Heyer M., Rathborne J., Simon R., 2010, ApJ, 723, 492
  * Sandstrom et al. (2013) Sandstrom K. M., Leroy A. K., Walter F., Bolatto A. D., Croxall K. V., Draine B. T., Wilson C. D., Wolfire M., Calzetti D., Kennicutt R. C., Aniano G., et al. 2013, ApJ, 777, 5
  * Schöier et al. (2005) Schöier F. L., van der Tak F. F. S., van Dishoeck E. F., Black J. H., 2005, A&A, 432, 369
  * Sembach et al. (2000) Sembach K. R., Howk J. C., Ryans R. S. I., Keenan F. P., 2000, ApJ, 528, 310
  * Shetty et al. (2012) Shetty R., Beaumont C. N., Burton M. G., Kelly B. C., Klessen R. S., 2012, MNRAS, 425, 720
  * Shetty et al. (011a) Shetty R., Glover S. C., Dullemond C. P., Klessen R. S., 2011a, MNRAS, 412, 1686
  * Shetty et al. (011b) Shetty R., Glover S. C., Dullemond C. P., Ostriker E. C., Harris A. I., Klessen R. S., 2011b, MNRAS, 415, 3253
  * Smith et al. (014a) Smith R. J., Glover S. C. O., Clark P. C., Klessen R. S., Springel V., 2014a, MNRAS, 441, 1628
  * Smith et al. (014b) Smith R. J., Glover S. C. O., Klessen R. S., 2014b, MNRAS, 445, 2900
  * Sobolev (1957) Sobolev V. V., 1957, SvA, 1, 678
  * Solomon et al. (1987) Solomon P. M., Rivolo A. R., Barrett J., Yahil A., 1987, ApJ, 319, 730
  * Springel (2010) Springel V., 2010, MNRAS, 401, 791
  * Szűcs et al. (2014) Szűcs L., Glover S. C. O., Klessen R. S., 2014, MNRAS, 445, 4055
  * Truelove et al. (1998) Truelove J. K., Klein R. I., McKee C. F., Holliman II J. H., Howell L. H., Greenough J. A., Woods D. T., 1998, ApJ, 495, 821
  * van der Tak & van Dishoeck (2000) van der Tak F. F. S., van Dishoeck E. F., 2000, A&A, 358, L79
  * Visser et al. (2009) Visser R., van Dishoeck E. F., Black J. H., 2009, A&A, 503, 323
  * Watson et al. (1976) Watson W. D., Anicich V. G., Huntress Jr. W. T., 1976, ApJ, 205, L165
  * Young & Scoville (1991) Young J. S., Scoville N. Z., 1991, ARA&A, 29, 581
  * Yusef-Zadeh et al. (2007) Yusef-Zadeh F., Wardle M., Roy S., 2007, ApJ, 665, L123

##  Appendix A Radiative transfer post-processing: influence of grid
resolution

Figure 9: Emission PDFs for our various tracers for different resolution
models with $128^{3}$, $256^{3}$ and $512^{3}$ grid cells for clouds with
$\alpha=2.0$.

Since Radmc-3d cannot deal with Arepo data directly, we have to map the
simulation output onto a cubic grid. The results presented in this paper are
all based on a grid resolution of 5123 cells. However, we examine the
sensitivity of our results to the choice of the number of grid cells in each
dimension. Therefore, we run a similar radiative transfer post-processing with
Radmc-3d, using resolutions of 2563 and 1283 grid zones. In Fig. 9 we compare
the probability density functions (PDF) of the velocity-integrated intensities
for our various tracers for the different grid resolutions, while using the
intermediate $\alpha=2.0$ cloud models as an example. We find that a
resolution of $512^{3}$ grid cells is enough in order to properly recover the
emission of the cloud for all of the tracers apart from the [Oi] 63$\mu$m
line. In the case of [Oi] 63$\mu$m, the maps are well-converged for integrated
intensities of below $1.5\,{\rm K\,km\,s^{-1}}$, but at higher integrated
intensities, increasing the resolution depresses the PDF. However, as we have
already seen in Section 3.4, integrated intensities of this magnitude are only
recovered for [Oi] 63$\mu$m along lines of sight where the line is already
optically thick. The precise values of the integrated intensity that we
recover in this regime therefore do not significantly affect our conclusions.

For all of the tracers, we find that the differences between the PDFs we
recover from the $256^{3}$ and $512^{3}$ runs are always significantly smaller
than the differences between the $128^{3}$ and $256^{3}$ models. Thus, we do
not expect to find larger differences in the PDFs for models with even higher
resolution, e.g. between $512^{3}$ and $1024^{3}$. Hence, due to the
significant increase of the computational costs for the radiative transfer
post-processing for even higher resolution runs, we focus on models with a
resolution of $512^{3}$ grid cells in this paper, which is enough to
demonstrate our basic conclusions.

##  Appendix B Dependence on viewing angle

Figure 10: Example of emission PDFs for the 12CO tracer for the resolution
model with $512^{3}$ grid cells. We choose the intermediate cloud model with
$\alpha=2.0$. The above figure shows the emission PDF for the $x$-, $y$\- and
$z$-direction. Although there are slight variations in the different PDFs, the
bulk emission is the same for all directions.

When constructing our synthetic emission maps, we choose to focus on a line of
sight (LoS) parallel to the $z$-axis of the simulation volume. However, as the
turbulence in our simulations is isotropic, we expect our results to be
insensitive to this choice. To test this, we have also made maps of 12CO for
LoS parallel to the $x$ and $y$ axes, using our standard grid resolution of
5123 zones, for the run with $\alpha=2.0$. The resulting integrated intensity
PDFs are shown in Figure 10. Although there are slight variations in the
different PDFs for $x$, $y$ and $z$, the bulk emission is the same for all
directions. We find a similar behavior also for the other atomic tracers.
Hence, we conclude that it is enough to focus on one specific LoS, for example
the $z$-direction, in order to get an idea about the underlying physical
parameters of the cloud. We note, however, that for magnetized clouds this is
not necessarily true, as in this case the turbulence will no longer be
isotropic, unless the field is very weak.

##  Appendix C Influence of the CO abundance on the line emission of [Cii] and
[Oi]

Figure 11: Emission PDFs for our various atomic tracers for a fixed resolution
of $512^{3}$ grid cells for clouds with $\alpha=2.0$. Shown are the PDFs given
in Fig. 9 (with CO) together with PDFs computed from additional runs of
radiative transfer post-processing with updated abundances of carbon and
oxygen assuming zero CO abundances in our simulations (without CO). We do not
find significant changes in the different emission PDFs of the atomic tracers.

In Section 2.1 we emphasize that the simplified Nelson & Langer (1997) network
that we use in our study is known to somewhat overestimate the rate at which
CO forms in our simulations. Thus, we may expect that this has a consequence
on the distribution of the line emission of [Cii] and [Oi]. In order to
quantify this, we perform additional runs of radiative transfer post-
processing with updated abundances of carbon and oxygen for comparison by
assuming zero CO abundances in our simulations. As an example, we plot the
different line emission PDFs for our various atomic tracers for a fixed
resolution of $512^{3}$ grid cells for clouds with $\alpha=2.0$. The result is
shown in Fig. 11. We find that the [Cii] 158 $\mu$m and [Oi] 63 $\mu$m line
emission are unaffected by the changed abundances. This suggests that in the
region where the CO abundance is uncertain, these two lines are already
optically thick and hence are in sensitive to this uncertainty. Support for
this interpretation comes from the behavior of the [Oi] 145 $\mu$m line. As we
have already seen, this has a much smaller optical depth, and hence is more
strongly altered by the changed abundances. However, even in this case, the
effect is relatively small: the brightest regions become around 0.1 dex
brighter, but the majority of the PDF remains unaffected. We can therefore
conclude that the inaccuracy introduced into our predicted [Cii] and [Oi] maps
by our simplified chemical model is unimportant.


# Viscous dynamics of vortices in a ferromagnetic film

Derek Reitz  Anirban Ghosh  Oleg Tchernyshyov  Department of Physics and
Astronomy, Johns Hopkins University, Baltimore, Maryland 21218, USA

###### Abstract

We derive viscous forces for vortices in a thin-film ferromagnet. The viscous
force acting on vortex $i$ is a linear superposition
$\mathbf{F}_{i}=-\sum_{j}\hat{D}_{ij}\mathbf{V}_{j}$, where $\mathbf{V}_{j}$
is the velocity of vortex $j$. Thanks to the long-range nature of vortices,
the mutual drag tensor $\hat{D}_{ij}$ is comparable in magnitude to the
coefficient of self-drag $D_{ii}$.

The dynamics of solitons in ferromagnets is a topic with a long history. Time
evolution of magnetization, represented by the field
$\mathbf{m}(\mathbf{r},t)$ of unit length, is described by the Landau-Lifshitz
equation

$\mathcal{J}\dot{\mathbf{m}}=-\mathbf{m}\times\frac{\delta
U}{\delta\mathbf{m}}-\alpha|\mathcal{J}|\mathbf{m}\times\dot{\mathbf{m}},$ (1)

where $U[\mathbf{m}(\mathbf{r})]$ is a functional of potential energy,
$\mathcal{J}$ is the density of angular momentum sig , and $\alpha\ll 1$ is
Gilbert’s damping constant Gilbert (2004). Even in the simplest models, where
the energy includes only exchange interactions and local anisotropy, Eq. (1)
is a nonlinear partial differential equation that rarely admits exact
dynamical solutions. Approximate solutions can be obtained for soft modes
associated with global symmetries (such as translations) in the limit of weak
external perturbations. Thiele (1973) described the dynamics of a rigidly
moving magnetic soliton, $\mathbf{m}(\mathbf{r}-\mathbf{R}(t))$, whose
velocity $\dot{\mathbf{R}}$ is determined from the equation

$\mathbf{G}\times\dot{\mathbf{R}}-\partial
U/\partial\mathbf{R}-D\dot{\mathbf{R}}=0,$ (2)

expressing the balance of gyroscopic, potential, and viscous forces,
respectively.

Figure 1: (a) Vortex–antivortex pair with equal skyrmion numbers. Red and blue
colors signify positive and negative out-of-plane magnetization $m_{z}$,
respectively. (b) Trajectory of the vortex. Natural units of length.

Thiele’s equation (2) has been widely used to describe the dynamics of
vortices in a thin film Takeno and Homma (1982); Huber (1982); Nikiforov and
Sonin (1983); Mertens  _et al._ (1997); Mertens and Bishop (2000); Gaididei
_et al._ (2000); Guslienko  _et al._ (2002); Kovalev  _et al._ (2002);
Buchanan  _et al._ (2005); Sheka  _et al._ (2005); Komineas (2007); Ivanov
_et al._ (2010). With few exceptions Huber (1982); Mertens and Bishop (2000);
Sheka  _et al._ (2005), analytical treatments take into account the gyroscopic
and potential forces but leave out the viscous force involved in energy
dissipation. Both the gyroscopic and viscous forces are proportional to the
soliton velocity $\dot{\mathbf{R}}$ and the neglect of the viscous force can
be justified by its relative weakness: the viscosity tensor
$D_{ab}=\alpha|\mathcal{J}|\int
dV\,\partial_{a}\mathbf{m}\cdot\partial_{b}\mathbf{m}$ is of a higher order in
$\alpha\ll 1$ than the gyrovector $G_{a}=\epsilon_{abc}\mathcal{J}\int
dV\,\mathbf{m}\cdot(\partial_{b}\mathbf{m}\times\partial_{c}\mathbf{m})$.
Nonetheless, in certain situations the viscous force cannot be neglected. For
example, the annihilation of a vortex and an antivortex is accompanied by
gradual dissipation of energy as the two solitons approach each other. This
motivates us to seek a proper understanding of viscous forces in vortex
dynamics.

To be specific, we set as our immediate goal to obtain a satisfactory
analytical model for the motion of a vortex–antivortex pair with equal
skyrmion numbers, Fig. 1(a). The two solitons attract each other through a
potential force mediated by exchange interaction. To the zeroth order in
$\alpha$, Thiele’s equation (2), applied to each soliton separately, predicts
that they will orbit a common center at an orbital velocity proportional to
the force of attraction. At the next order in $\alpha$, viscous forces
opposing the orbital motion will induce slow radial motion of the solitons
toward each other, Fig. 1(b).

Our main findings are as follows. (1) Viscous forces acting on vortices come
in two kinds. The first is self-drag, a force proportional to the vortex’s own
velocity Huber (1982); Mertens and Bishop (2000). We show that a vortex also
experiences a drag force from other vortices proportional to _their_
velocities. The net force of viscous friction acting on vortex $i$ is
$\mathbf{F}_{i}=-\sum_{j}\hat{D}_{ij}\mathbf{V}_{j}$, where $\mathbf{V}_{j}$
is the velocity of vortex $j$ and $\hat{D}_{ij}$ is the mutual drag tensor
comparable in magnitude to the self-drag coefficient $D_{ii}$. (2) Both
$D_{ii}$ and $\hat{D}_{ij}$ scale logarithmically with the system size. (3)
The direction of mutual drag depends on the product of vorticities. A vortex
receding from an antivortex attempts to drag the antivortex with it; the force
direction reverses for a vortex–vortex pair. (4) Image vortices, created by
“reflection” in the sample edge, produce substantial corrections to viscous
forces.

Our theory is built on the framework of collective coordinates Mertens and
Bishop (2000); Tretiakov  _et al._ (2008); Clarke  _et al._ (2008), in which
the magnetization field $\mathbf{m}(\mathbf{r},t)$ is parametrized by a few
coordinates $\\{q^{\mu}\\}$, $\mu=1,2,\ldots$, representing soft modes of the
system (e.g., vortex positions). The Landau-Lifshitz equation (1) translates
into equations of motion for each coordinate $q^{\mu}$,

$G_{\mu\nu}\dot{q}^{\nu}-\frac{\partial U}{\partial
q^{\mu}}-D_{\mu\nu}\dot{q}^{\nu}=0.$ (3)

Thiele’s equation (2) is a particular case of Eq. (3), in which
$\\{q^{\mu}\\}$ are global translations
$\mathbf{m}(\mathbf{r})\mapsto\mathbf{m}(\mathbf{r}-\mathbf{R})$. The
gyroscopic and dissipative tensors are Tretiakov  _et al._ (2008); Clarke  _et
al._ (2008)

$\displaystyle G_{\mu\nu}$ $\displaystyle=$ $\displaystyle-\mathcal{J}\int
dV\,\mathbf{m}\cdot\left(\frac{\partial\mathbf{m}}{\partial
q^{\mu}}\times\frac{\partial\mathbf{m}}{\partial q^{\nu}}\right),$ (4a)
$\displaystyle D_{\mu\nu}$ $\displaystyle=$
$\displaystyle\alpha|\mathcal{J}|\int dV\,\frac{\partial\mathbf{m}}{\partial
q^{\mu}}\cdot\frac{\partial\mathbf{m}}{\partial q^{\nu}}.$ (4b)

We use a simple model of a thin-film ferromagnet with exchange interaction of
strength $A$ and easy-plane anisotropy of strength $K$. We omit dipolar
interactions dip . In a film of thickness $h$, the energy is

$U=h\int d^{2}r\left(A|\nabla\mathbf{m}|^{2}+Km_{z}^{2}\right)/2.$ (5)

The polar angle of magnetization $\theta$ is a hard mode pinned at
$\theta=\pi/2$. At low energies, the system is effectively an $XY$ ferromagnet
Chaikin and Lubensky (2000) parametrized by the azimuthal angle of
magnetization $\phi(\mathbf{r},t)$ with energy $U=h\int
d^{2}r\,A|\nabla\phi|^{2}/2$. A state with $N$ vortices has

$\phi(\mathbf{r})=\sum_{i=1}^{N}n_{i}\arctan{\frac{y-Y_{i}}{x-X_{i}}},$ (6)

where $\mathbf{R}_{i}=(X_{i},Y_{i})$ is the location of the $i$th vortex and
$n_{i}\in\mathbb{Z}$ is its vorticity, usually $n_{i}=\pm 1$. The effective
description breaks down inside vortex cores—circular regions with the size on
the scale of the exchange length $\lambda=\sqrt{A/K}$, where $\mathbf{m}$
comes out of the easy plane.

Gyroscopic and potential forces acting on vortices are well understood. The
gyroscopic density in (4a) comes from core regions, where $\mathbf{m}$ does
not stay in a fixed plane Clarke  _et al._ (2008). Vortex cores are rigid
objects, for which Thiele’s approximation Thiele (1973) works well. The
gyroscopic force $\mathbf{F}^{g}$ for a vortex with velocity
$\dot{\mathbf{R}}=(\dot{X},\dot{Y})$ has components

$F^{g}_{X}=-4\pi Q\mathcal{J}h\dot{Y},\quad F^{g}_{Y}=4\pi
Q\mathcal{J}h\dot{X},$ (7)

where $Q=np/2=\pm 1/2$ is the skyrmion number of the vortex determined by its
vorticity $n$ and polarity $p=\pm 1$. Exchange-mediated conservative forces
between vortices resemble Coulomb interactions in two dimensions Kosterlitz
(1974). The net conservative force on vortex $i$ is

$\mathbf{F}^{c}_{i}=-2\pi Ah\sum_{j\neq
i}n_{i}n_{j}\frac{\mathbf{R}_{i}-\mathbf{R}_{j}}{|\mathbf{R}_{i}-\mathbf{R}_{j}|^{2}}.$
(8)

Viscous forces are the primary focus of this paper. It is natural to expect
that vortex $i$ experiences a viscous force
$\mathbf{F}^{v}_{i}=-\hat{D}_{ii}\dot{\mathbf{R}}_{i}$, where $\hat{D}_{ii}$
is a $2\times 2$ symmetric tensor with matrix elements such as
$D_{X_{i}X_{i}}$, $D_{X_{i}Y_{i}}$, and so on. For magnetization lying
primarily in the easy plane, $\mathbf{m}\approx(\cos{\phi},\sin{\phi},0)$, Eq.
(4b) yields, e.g.,

$D_{X_{i}X_{i}}\approx\alpha|\mathcal{J}|h\int
d^{2}r\left(\frac{\partial\phi}{\partial
X_{i}}\right)^{2}=\alpha|\mathcal{J}|h\int
d^{2}r\frac{(y-Y_{i})^{2}}{|\mathbf{r}-\mathbf{R}_{i}|^{4}}.$

On symmetry grounds, we expect $\hat{D}_{ii}$ to be isotropic,
$\hat{D}_{ii}=D_{ii}\hat{1}$, with a scalar viscosity coefficient

$D_{ii}=\frac{\alpha|\mathcal{J}|h}{2}\int\frac{d^{2}r}{|\mathbf{r}-\mathbf{R}_{i}|^{2}}.$
(9)

The integral diverges and requires regularization for both
$\mathbf{r}\to\mathbf{R}_{i}$ and $\mathbf{r}\to\infty$. The long-range cutoff
is the system size Huber (1982); the short-range cutoff is provided by the
size of the vortex core of the order $\lambda$ Sheka  _et al._ (2005). For a
vortex near the center of a disk of radius $R_{d}\gg R_{i}$,

$D_{ii}\approx\alpha\pi|\mathcal{J}|h\ln{(R_{d}/C\lambda)},$ (10)

where $C$ is a numerical factor of the order 1. The logarithmic divergence
with the system size $R_{d}$ reflects the long-range impact of a moving vortex
on the magnetization distribution $\mathbf{m}(\mathbf{r})$. Viscous forces in
a ferromagnet are of the order $\alpha\ll 1$ (typically $10^{-4}$ to
$10^{-2}$) and thus are much weaker than gyroscopic ones. For a vortex, this
is partly compensated by the factor $\ln{(R_{d}/C\lambda)}\gg 1$.

The extended nature of vortices leads to substantial _mutual drag_ between
them. Vortex $i$ feels a force proportional to the velocity of vortex $j$,
$\mathbf{F}^{v}_{i}=-\hat{D}_{ij}\dot{\mathbf{R}}_{j}$, where again
$\hat{D}_{ij}$ is a $2\times 2$ tensor with coefficients such as

$\displaystyle D_{X_{i}X_{j}}$ $\displaystyle\approx$
$\displaystyle\alpha|\mathcal{J}|h\int d^{2}r\,\frac{\partial\phi}{\partial
X_{i}}\frac{\partial\phi}{\partial X_{j}}$ $\displaystyle=$
$\displaystyle\alpha|\mathcal{J}|hn_{i}n_{j}\int
d^{2}r\frac{(y-Y_{i})(y-Y_{j})}{|\mathbf{r}-\mathbf{R}_{i}|^{2}|\mathbf{r}-\mathbf{R}_{j}|^{2}}.$

The integrand has two singularities at $\mathbf{r}=\mathbf{R}_{i}$ and
$\mathbf{R}_{j}$. They are weaker than the confluent singularity in Eq. (9)
and are integrable, making the short-range cut-off $\lambda$ unnecessary. For
two vortices located symmetrically about the center of a disk,
$\mathbf{R}_{i}=(\pm R,0)$ in Fig. 1, the viscosity tensor $\hat{D}_{ij}$ has
principal axes parallel and perpendicular to the line connecting the cores
with the eigenvalues

$\displaystyle D_{ij}^{||}$ $\displaystyle\approx$ $\displaystyle\alpha\pi
n_{i}n_{j}|\mathcal{J}|h\left[\ln{(R_{d}/2R)}+1/2\right],$ (11a)
$\displaystyle D_{ij}^{\perp}$ $\displaystyle\approx$ $\displaystyle\alpha\pi
n_{i}n_{j}|\mathcal{J}|h\left[\ln{(R_{d}/2R)}-1/2\right].$ (11b)

Table 1: Vortices in the numerical simulation. $i$ | $n_{i}$ | $Q_{i}$ | $G_{X_{i}Y_{i}}$ | $R_{i}$ | $\Phi_{i}$ | note
---|---|---|---|---|---|---
1 | $+1$ | $+1/2$ | $-2\pi\mathcal{J}$ | $R$ | $\Phi$ | vortex
2 | $-1$ | $+1/2$ | $-2\pi\mathcal{J}$ | $R$ | $\Phi+\pi$ | antivortex
3 | $-1$ |  | 0 | $R_{d}^{2}/R$ | $\Phi$ | image of the vortex
4 | $+1$ |  | 0 | $R_{d}^{2}/R$ | $\Phi+\pi$ | image of the antivortex

A notable feature of mutual drag is the dependence of its direction on
vorticities. The drag force
$\mathbf{F}^{v}_{i}=-\hat{D}_{ij}\dot{\mathbf{R}}_{j}$ on vortex $i$ is
roughly opposite to the velocity $\dot{\mathbf{R}}_{j}$ of vortex $j$ for
vorticities of the same sign and roughly parallel to it for vorticities of
opposite sign.

To test our theory, we modeled the dynamics of a vortex–antivortex pair with
equal skyrmion numbers (Table 1) with the aid of the micromagnetic simulator
MuMax3 Vansteenkiste  _et al._ (2014). We used magnetization length
$\mathcal{M}=8.60\times 10^{5}\ \mbox{A/m}$, gyromagnetic ratio
$\gamma=-2.21\times 10^{5}\ \mbox{m/A s}$, angular momentum density
$\mathcal{J}=\mu_{0}\mathcal{M}/\gamma=-4.89\times 10^{-6}$ J s/m3 sig ,
exchange constant $A=2.6\times 10^{-11}\ \mbox{J/m}$, and easy-plane
anisotropy $K=2.60\times 10^{5}\ \mbox{J/m}^{3}$. Natural units of length and
time were $\lambda=\sqrt{A/K}=10.0\mbox{ nm}$ and $\tau=|\mathcal{J}|/K=18.8$
ps. The sample was a disk of radius $R_{D}=2048$ nm and thickness $h=4$ nm. A
vortex and an antivortex with equal skyrmion numbers $Q_{1}=Q_{2}=+1/2$ were
initially placed symmetrically on opposite sides of the disk center, Fig.
1(a). The pair orbited the disk center and gradually spiraled down, Fig. 1(b).
The dimensionless constant $C=0.342$ in Eq. (10) was determined through a
numerical evaluation of the dissipation constant
$D_{ii}=\alpha|\mathcal{J}|h\int d^{2}r|\partial_{x}\mathbf{m}|^{2}$ of a
simulated vortex.

In the absence of dissipation ($\alpha=0$), the motion of the pair reflects
the balance of the exchange-mediated attraction (8) and the gyroscopic force
(7), $\mathbf{F}^{c}+\mathbf{F}^{g}=0$. The radial direction of the exchange
attraction results in the azimuthal direction of the vortex velocities. The
two topological defects orbit the common center at a constant radius $R$. It
is therefore convenient to parametrize the positions of the vortices in polar
coordinates $(R_{i},\Phi_{i})$, see Table 1. The angular velocity is obtained
from the balance of the gyroscopic and conservative forces acting on a vortex
in the radial direction,

$-2\pi\mathcal{J}hR\dot{\Phi}-\frac{2\pi Ah}{2R}=0.$ (12)

Weak dissipation ($\alpha\ll 1$) turns the trajectories into spirals with a
radial velocity $\dot{R}$ of the order $\alpha$.

Figure 2: (a) Angular velocity $\dot{\Phi}$ of the vortices vs. the radius $R$
of their orbit. (b) The ratio of radial $\dot{R}$ and orbital $R\dot{\Phi}$
velocities vs. the orbit radius $R$. Natural units of length $\lambda$ and
time $\tau$.

Numerical simulations reveal very good, but not perfect, agreement with Eq.
(12): the observed angular velocity of the vortices $\dot{\Phi}(R)$ differed
from the expected value $\dot{\Phi}=-A/2\mathcal{J}R^{2}$ by a small constant,
Fig. 2(a). This minor discrepancy reflects an edge effect in a finite system.
Free boundary conditions at the edge, $\partial\phi/\partial r=0$, result in
the appearance of image vortices outside of the disk Gaididei  _et al._
(2000); Sheka  _et al._ (2005), see Table 1. The images generate a weak radial
force of approximately $4\pi AhR/R_{d}^{2}$ that reduces the angular velocity
by $2A/\mathcal{J}R_{d}^{2}$, in excellent agreement with the numerical data,
Fig. 2(a).

To determine the radial velocity $\dot{R}$, which is of the order $\alpha$, we
need to carefully evaluate viscous forces acting on the vortices. This task is
made complicated by the constrained motion of images (their positions mirror
the locations of the vortices). Because of these constraints, the vortices
also “feel” forces acting on the images. Although it is possible to solve the
dynamics with constraints, a more expedient way is to reformulate the dynamics
in terms of the two independent variables $R$ and $\Phi$ that fully determine
the positions of all four objects (see Table 1). Equations of motion for $R$
and $\Phi$ can be obtained by following the usual prescription (3):

$\displaystyle G_{R\Phi}\dot{\Phi}-\partial U/\partial R-D_{RR}\dot{R}=0,$
(13a) $\displaystyle G_{\Phi R}\dot{R}-D_{\Phi\Phi}\dot{\Phi}=0.$ (13b)

Here we took into account rotational symmetry, which yields $-\partial
U/\partial\Phi=0$, and the diagonal nature of the dissipation tensor,
$D_{R\Phi}=D_{\Phi R}=0$.

Polar components of the gyroscopic tensor $G_{R\Phi}=-G_{\Phi R}$ can be
expressed in terms of Cartesian ones through a standard coordinate change from
$\\{q^{\mu}\\}$ (here positions of the vortices and images $\mathbf{R}_{i}$)
to $R$ and $\Phi$:

$\displaystyle G_{R\Phi}$ $\displaystyle=$
$\displaystyle\sum_{\mu,\nu}\frac{\partial q^{\mu}}{\partial R}\frac{\partial
q^{\nu}}{\partial\Phi}G_{q^{\mu}q^{\nu}}=R(G_{X_{1}Y_{1}}+G_{X_{2}Y_{2}})$
(14) $\displaystyle=$ $\displaystyle-4\pi R\mathcal{J}h.$

Here we used the reference frame of Fig. 1, in which
$\partial\mathbf{R}_{1}/\partial R=(1,0)$,
$\partial\mathbf{R}_{1}/\partial\Phi=(0,R)$, $\partial\mathbf{R}_{2}/\partial
R=(-1,0)$, $\partial\mathbf{R}_{2}/\partial\Phi=(0,-R)$, etc. Image vortices 3
and 4 do not contribute to the gyroscopic tensor because they are centered
outside the sample and therefore lack cores. The energy of the vortices and
images is

$U(R)=2\pi
Ah\left(\ln{2R}+\ln{\frac{R_{d}^{2}-R^{2}}{R_{d}^{2}+R^{2}}}\right).$ (15)

The dissipative term in Eq. (13a) is of the order $\alpha^{2}$ and can be
neglected to yield $\dot{\Phi}\approx-(A/2\mathcal{J})(R^{-2}-4R_{d}^{-2})$
for $R\ll R_{d}$, as derived above.

To obtain the radial velocity from Eq. (13b), we need the dissipative
coefficient

$D_{\Phi\Phi}=\sum_{\mu,\nu}\frac{\partial
q^{\mu}}{\partial\Phi}\frac{\partial
q^{\nu}}{\partial\Phi}D_{q^{\mu}q^{\nu}}=\sum_{i,j}(-1)^{i+j}R_{i}R_{j}D_{ij}^{\perp},$
(16)

which reduces to a superposition of these terms:

$\displaystyle 2R_{1}^{2}D_{11}\approx
2\alpha\pi|\mathcal{J}|hR^{2}\ln{(R_{d}/C\lambda)},$ (17a)
$\displaystyle-2R_{1}^{2}D_{12}^{\perp}\approx
2\alpha\pi|\mathcal{J}|hR^{2}[\ln{(R_{d}/2R)}-1/2],$ (17b) $\displaystyle
4R_{1}R_{3}(D_{13}^{\perp}-D_{14}^{\perp})\approx
4\alpha\pi|\mathcal{J}|hR^{2},$ (17c) $\displaystyle
2R_{3}^{2}(D_{33}^{\perp}-D_{34}^{\perp})\approx\alpha\pi|\mathcal{J}|hR^{2}.$
(17d)

They represent self-drag of the vortices (17a), their mutual drag (17b), and
corrections from the images (17c) and (17d). In the limit of a large disk, the
first two terms dominate over the edge corrections, albeit only
logarithmically in the system size $R_{d}$. Thus the edge corrections must be
included to obtain quantitative agreement with simulations, Fig. 2(b).

The ratio of the radial and orbital velocities,

$\frac{\dot{R}}{R\dot{\Phi}}=\frac{D_{\Phi\Phi}}{RG_{\Phi
R}}=-\mathrm{sgn}{\mathcal{J}}\frac{\alpha}{2}\ln{\frac{\Lambda}{R}},$ (18)

where $\Lambda=R_{d}^{2}e^{2}/2C\lambda$, has a telltale logarithmic
dependence on the vortex separation $2R$ inherited from mutual drag, Eqs.
(11b) and (17b). As a result, a vortex-antivortex pair follows a _double_
logarithmic spiral,

$\ln{\ln{\frac{R(\Phi)}{\Lambda}}}-\ln{\ln{\frac{R(0)}{\Lambda}}}=-\mathrm{sgn}{\mathcal{J}}\frac{\alpha\Phi}{2},$
(19)

in excellent agreement with the numerical simulations, Fig. 1(b). Note the
contrast with a single vortex in a disk, which has a constant self-drag
coefficient (10) and therefore follows a simple logarithmic spiral Sheka  _et
al._ (2005).

At the smallest orbital radii $R$, the observed radial velocity $\dot{R}$
shows small but growing deviations from the theoretical value, Fig. 2(b). As
the radial motion is tied to energy dissipation, an excess radial velocity
hints at the opening of a new dissipation channel. The likely culprit is spin
waves, which have a linear dispersion $\omega=sk$ with the speed
$s=\sqrt{AK/\mathcal{J}^{2}}=\lambda/\tau$. In a disk, the normal modes in
polar coordinates $(R,\Phi)$ are

$\phi(R,\Phi,t)=aJ_{m}(kR)\cos{(\omega t-m\Phi)},$ (20)

where $J_{m}(x)$ is a Bessel function of the first kind. For open boundary
conditions, the wavenumbers $k$ satisfy $J_{m}^{\prime}(kR_{d})=0$. A rotating
vortex-antivortex pair couples strongly to modes with $m=1$. The lowest
frequency for an $m=1$ spin wave is $\omega=1.84s/R_{d}=8.98\times
10^{-3}\tau^{-1}$. At the end of the simulation, the angular frequency of the
pair reached $\dot{\Phi}=4.3\times 10^{-3}\tau^{-1}$. Although the pair was
not yet in resonance with this mode, its angular velocity had a substantial
chirp, $|\ddot{\Phi}/\dot{\Phi}^{2}|=\alpha\ln{(\Lambda/R)}\approx 0.1$, and
thus a spectrum potentially wide enough to excite the $m=1$ spin wave, whose
dynamics would produce additional dissipation.

We have derived viscous forces acting on a vortex in a thin-film ferromagnet.
In addition to self-drag proportional to the vortex’s own velocity, vortices
experience mutual drag, a force on vortex $i$ proportional to the velocity of
vortex $j$, $\mathbf{F}^{v}_{i}=-\hat{D}_{ij}\dot{\mathbf{R}}_{j}$. Reflecting
the long-range influence of vortices, both the self-drag coefficient $D_{ii}$
(10) and the mutual viscosity tensor $\hat{D}_{ij}$ (11) scale logarithmically
with the system size. The mutual drag tensor $\hat{D}_{ij}$ is anisotropic and
distance-dependent. We have tested our theory by deriving the dynamics of a
vortex-antivortex pair with equal skyrmion numbers. The results are in
excellent agreement with micromagnetic simulations. Edge effects in the form
of image vortices contribute substantially to viscous friction. Dissipation
through the emission of spin waves becomes noticeable when vortices approach
each other very closely, within a few exchange lengths $\lambda$.

The authors thank Se Kwon Kim for helpful discussions. This work was supported
by the US Department of Energy, Office of Basic Energy Sciences, Division of
Materials Sciences and Engineering under Award DE-FG02-08ER46544.

## References

  * (1) Here $\mathcal{J}=\mu_{0}\mathcal{M}/\gamma$, where $\mathcal{M}>0$ is magnetization length and $\gamma$ is the gyromagnetic ratio. Thus $\mathrm{sgn}\mathcal{J}=\mathrm{sgn}\gamma$.
  * Gilbert (2004) T. L. Gilbert, IEEE Trans. Mag. 40, 3443 (2004).
  * Thiele (1973) A. A. Thiele, Phys. Rev. Lett. 30, 230 (1973).
  * Takeno and Homma (1982) S. Takeno and S. Homma, Prog. Theor. Phys. 67, 1633 (1982).
  * Huber (1982) D. L. Huber, Phys. Rev. B 26, 3758 (1982).
  * Nikiforov and Sonin (1983) A. V. Nikiforov and É. B. Sonin, Sov. Phys. JETP 58, 373 (1983).
  * Mertens  _et al._ (1997) F. G. Mertens, H. J. Schnitzer, and A. R. Bishop, Phys. Rev. B 56, 2510 (1997).
  * Mertens and Bishop (2000) F. G. Mertens and A. R. Bishop, in  _Nonlinear Science at the Dawn of the 21th Century_, Lecture Notes in Physics, edited by P. L. Christiansen, M. P. Soerensen, and A. C. Scott (Springer, Berlin, 2000) pp. 137–170.
  * Gaididei  _et al._ (2000) Y. Gaididei, T. Kamppeter, F. G. Mertens, and A. R. Bishop, Phys. Rev. B 61, 9449 (2000).
  * Guslienko  _et al._ (2002) K. Y. Guslienko, B. A. Ivanov, V. Novosad, Y. Otani, H. Shima, and K. Fukamichi, J. Appl. Phys. 91, 8037 (2002).
  * Kovalev  _et al._ (2002) A. Kovalev, S. Komineas, and F. Mertens, Eur. J. Phys. B 25, 89 (2002).
  * Buchanan  _et al._ (2005) K. S. Buchanan, P. E. Roy, M. Grimsditch, F. Y. Fradin, K. Y. Guslienko, S. D. Bader, and V. Novosad, Nat. Phys. 1, 172 (2005).
  * Sheka  _et al._ (2005) D. D. Sheka, J. P. Zagorodny, J.-G. Caputo, Y. Gaididei, and F. G. Mertens, Phys. Rev. B 71, 134420 (2005).
  * Komineas (2007) S. Komineas, Phys. Rev. Lett. 99, 117202 (2007).
  * Ivanov  _et al._ (2010) B. A. Ivanov, E. G. Galkina, and A. Y. Galkin, Low Temp. Phys. 36, 747 (2010).
  * Tretiakov  _et al._ (2008) O. A. Tretiakov, D. Clarke, G.-W. Chern, Y. B. Bazaliy, and O. Tchernyshyov, Phys. Rev. Lett. 100, 127204 (2008).
  * Clarke  _et al._ (2008) D. J. Clarke, O. A. Tretiakov, G.-W. Chern, Y. B. Bazaliy, and O. Tchernyshyov, Phys. Rev. B 78, 134412 (2008).
  * (18) To be exact, effects of dipolar interactions are partly included in our model. In a thin film, the anisotropy constant $K$ includes a dipolar contribution $\mu_{0}\mathcal{M}^{2}/2$. We neglect the nonlocal stray-field effects.
  * Chaikin and Lubensky (2000) P. M. Chaikin and T. C. Lubensky,  _Principles of Condensed Matter Physics_ (Cambridge University Press, 2000).
  * Kosterlitz (1974) J. M. Kosterlitz, J. Phys. C 7, 1046 (1974).
  * Vansteenkiste  _et al._ (2014) A. Vansteenkiste, J. Leliaert, M. Dvornik, M. Helsen, F. Garcia-Sanchez, and B. Van Waeyenberge, AIP Adv. 4, 107133 (2014).


* [39] M. Meyries and M.C. Veraar.  Sharp embedding results for spaces of smooth functions with power weights.  Studia Math., 208(3):257–293, 2012.
  * [40] M. Meyries and M.C. Veraar.  Traces and embeddings of anisotropic function spaces.  Math. Ann., 360(3-4):571–606, 2014.
  * [41] J.M.A.M. van Neerven, M.C. Veraar, and L.W. Weis.  Stochastic integration in UMD Banach spaces.  Ann. Probab., 35(4):1438–1478, 2007.
  * [42] J.M.A.M. van Neerven, M.C. Veraar, and L.W. Weis.  Stochastic maximal $L^{p}$-regularity.  Ann. Probab., 40(2):788–812, 2012.
  * [43] J. Prüss.  Evolutionary integral equations and applications.  Modern Birkhäuser Classics. Birkhäuser/Springer Basel AG, Basel, 1993.  [2012] reprint of the 1993 edition.
  * [44] J. Prüss.  Vector-valued Fourier multipliers in $L_{p}$-spaces with power weights.  Studia Math., 247(2):155–173, 2019.
  * [45] J. Prüss and G. Simonett.  Maximal regularity for evolution equations in weighted $L_{p}$-spaces.  Archiv der Mathematik, 82(5):415–431, 2004.
  * [46] J. Prüss and G. Simonett.  Moving interfaces and quasilinear parabolic evolution equations, volume 105.  Springer, 2016.
  * [47] J. Prüss, G. Simonett, and M. Wilke.  Critical spaces for quasilinear parabolic evolution equations and applications.  J. Differential Equations, 264(3):2028–2074, 2018.
  * [48] V.S. Rychkov.  On restrictions and extensions of the Besov and Triebel-Lizorkin spaces with respect to Lipschitz domains.  J. London Math. Soc. (2), 60(1):237–257, 1999.
  * [49] Y. Sawano.  Theory of Besov spaces, volume 56 of Developments in Mathematics.  Springer, Singapore, 2018.
  * [50] H.-J. Schmeißer and W. Sickel.  Traces, Gagliardo-Nirenberg Inequalities and Sobolev Type Embeddings for Vector-valued Function Spaces.  Unpublished notes, Jena, 2004.
  * [51] E. M. Stein.  Note on singular integrals.  Proceedings of the American Mathematical Society, 8(2):250–254, 1957.
  * [52] E.M. Stein.  Harmonic analysis: real-variable methods, orthogonality, and oscillatory integrals, volume 43 of Princeton Mathematical Series, Monographs in Harmonic Analysis, III.  Princeton University Press, Princeton, NJ, 1993.
  * [53] H. Triebel.  Interpolation theory, function spaces, differential operators, volume 18 of North-Holland Mathematical Library.  North-Holland Publishing Co., Amsterdam-New York, 1978.
  * [54] H. Triebel.  Theory of function spaces, volume 78 of Monographs in Mathematics.  Birkhäuser Verlag, Basel, 1983.
  * [55] H. Triebel.  Theory of function spaces. II, volume 84 of Monographs in Mathematics.  Birkhäuser Verlag, Basel, 1992.
  * [56] H. Triebel.  Fractals and spectra, volume 91 of Monographs in Mathematics.  Birkhäuser Verlag, Basel, 1997.  Related to Fourier analysis and function spaces.
  * [57] H. Triebel.  Theory of function spaces. III, volume 100 of Monographs in Mathematics.  Birkhäuser Verlag, Basel, 2006.
  * [58] R. Zacher.  Maximal regularity of type $L_{p}$ for abstract parabolic Volterra equations.  J. Evol. Equ., 5(1):79–103, 2005.