t_{1}}\otimes\frac{\partial\hat{\vec{j}}(\vec{x}_{2},t_{2})}{\partial
t_{2}}|0\rangle=\bm{1}\delta(x_{1}-x_{2})\int_{0}^{\infty}d\omega\int\frac{d^{2}\vec{k}_{\parallel}}{(2\pi)^{2}}\frac{\hbar\varepsilon_{0}}{\pi}\text{Im}[\varepsilon(\omega)]e^{i\vec{k}_{\parallel}\vec{\cdot}(\vec{x}_{1}-\vec{x}_{2})}\\\
\times\bigg{\\{}\Theta(x_{1}-a)\omega^{4}e^{-i\omega(t_{1}-t_{2})}+\Theta(-x_{1})\omega_{+}^{4}e^{-i\omega_{+}(t_{1}-t_{2})}\bigg{\\}}$
(154)

where $\omega_{\pm}=\omega\pm\vec{V}\vec{\cdot}\vec{k}_{\parallel}$, and
$\vec{k}_{\parallel}$ is a wave–vector lying in the $y$–$z$ plane. Inserting
the electric current correlation function (154) into (153) gives us the
electric field correlation function written without reference to the operators

$\langle
0|\hat{\vec{E}}(\vec{x},t)\otimes\hat{\vec{E}}(\vec{x}^{\prime},t)|0\rangle=\frac{\hbar\mu_{0}}{\pi
c^{2}}\int_{0}^{\infty}d\omega\int\frac{d^{2}\vec{k}_{\parallel}}{(2\pi)^{2}}\text{Im}[\varepsilon(\omega)]e^{i\vec{k}_{\parallel}\vec{\cdot}(\vec{x}-\vec{x}^{\prime})}\\\
\times\bigg{\\{}\omega^{4}\int_{a}^{\infty}dx_{1}\bm{G}(x,x_{1},\vec{k}_{\parallel},\omega)\vec{\cdot}\bm{G}^{\dagger}(x^{\prime},x_{1},\vec{k}_{\parallel},\omega)\\\
+\omega_{+}^{4}\int_{-\infty}^{0}dx_{1}\bm{G}(x,x_{1},\vec{k}_{\parallel},\omega_{+})\vec{\cdot}\bm{G}^{\dagger}(x^{\prime},x_{1},\vec{k}_{\parallel},\omega_{+})\bigg{\\}}$
(155)

where $\bm{G}^{\dagger}=(\bm{G}^{T})^{\star}$.

Integral identities for the Green function: To this order we are neglecting
all effects of the motion except the Doppler shift within the dispersion of
the medium. Therefore the differential equation satisfied by the Green
function is

$\vec{\nabla}\vec{\times}\vec{\nabla}\vec{\times}\bm{G}(x,x^{\prime},\vec{k}_{\parallel},\omega)-\frac{\omega^{2}}{c^{2}}[\varepsilon(\omega_{-})\Theta(-x)+\varepsilon(\omega)\Theta(x-a)]\bm{G}(x,x^{\prime},\vec{k}_{\parallel},\omega)=\bm{1}\delta(x-x^{\prime}),$
(156)

where $\vec{\nabla}=\vec{e}_{x}\partial_{x}+i\vec{k}_{\parallel}$, and
$\omega_{\pm}=\omega\pm\vec{V}\vec{\cdot}\vec{k}_{\parallel}$. Meanwhile the
Hermitian conjugate of the Green function obeys

$\vec{\nabla}\vec{\times}\vec{\nabla}\vec{\times}\bm{G}^{\dagger}(x^{\prime},x,\vec{k}_{\parallel},\omega)-\frac{\omega^{2}}{c^{2}}[\varepsilon^{\star}(\omega_{-})\Theta(-x)+\varepsilon^{\star}(\omega)\Theta(x-a)]\bm{G}^{\dagger}(x^{\prime},x,\vec{k}_{\parallel},\omega)=\bm{1}\delta(x-x^{\prime})$
(157)

Multiplying (157) on the left by the Green function,
$\bm{G}(x,x^{\prime\prime},\vec{k}_{\parallel},\omega)$, integrating over $x$,
and then subtracting the Hermitian conjugate of the resulting expression with
$x^{\prime}$ and $x^{\prime\prime}$ reversed, we obtain a generalisation of
(86),

$\frac{\bm{G}(x^{\prime},x^{\prime\prime},\vec{k}_{\parallel},\omega)-\bm{G}^{\dagger}(x^{\prime\prime},x^{\prime},\vec{k}_{\parallel},\omega)}{2i}=\frac{\omega^{2}}{c^{2}}\bigg{[}\text{Im}[\varepsilon(\omega_{-})]\int_{-\infty}^{0}dx\bm{G}(x^{\prime},x,\vec{k}_{\parallel},\omega)\vec{\cdot}\bm{G}^{\dagger}(x^{\prime\prime},x,\vec{k}_{\parallel},\omega)\\\
+\text{Im}[\varepsilon(\omega)]\int_{a}^{\infty}dx\bm{G}(x^{\prime},x,\vec{k}_{\parallel},\omega)\vec{\cdot}\bm{G}^{\dagger}(x^{\prime\prime},x,\vec{k}_{\parallel},\omega)\bigg{]}.$
(158)

For integration over only the region $x<0$ we have instead

$\frac{\omega^{2}}{c^{2}}\text{Im}[\varepsilon(\omega_{-})]\int_{-\infty}^{0}dx\bm{G}(x^{\prime},x,\vec{k}_{\parallel},\omega)\vec{\cdot}\bm{G}^{\dagger}(x^{\prime\prime},x,\vec{k}_{\parallel},\omega)=\frac{\bm{G}(x^{\prime},x^{\prime\prime},\vec{k}_{\parallel},\omega)-\bm{G}^{\dagger}(x^{\prime\prime},x^{\prime},\vec{k}_{\parallel},\omega)}{2i}\\\
-\frac{1}{2i}\left[\bm{G}(x^{\prime},0,\vec{k}_{\parallel},\omega)\vec{\cdot}\vec{e}_{x}\vec{\times}\vec{\nabla}\vec{\times}\bm{G}^{\dagger}(x^{\prime\prime},0,\vec{k}_{\parallel},\omega)-\text{h.c.}\right].$
(159)

Applying (158) to (155), gives us the electric field correlation function, and
after taking the limit, $\vec{x}\to\vec{x}^{\prime}$, we obtain the electric
contribution to the stress tensor

$\langle
0|\hat{\vec{E}}(\vec{x},t)\vec{\otimes}\hat{\vec{E}}(\vec{x},t)|0\rangle=\frac{\hbar\mu_{0}}{\pi}\int\frac{d^{2}\vec{k}_{\parallel}}{(2\pi)^{2}}\int_{0}^{\infty}d\omega\,\omega^{2}\text{Im}[\bar{\bm{G}}^{(S)}(x,x,\vec{k}_{\parallel},\omega)]\\\
-\int_{-\infty}^{0}dx_{1}\int\frac{d^{2}\vec{k}_{\parallel}}{(2\pi)^{2}}\int_{0}^{\vec{V}\vec{\cdot}\vec{k}_{\parallel}}d\omega\bar{\bm{Q}}(x,x_{1},\vec{k}_{\parallel},\omega),$
(160)

where we have introduced symmetrised quantities, e.g.
$\bar{\bm{G}}=\frac{1}{2}[\bm{G}+\bm{G}^{T}]$, which makes explicit the
independence of the right hand side from the order of the electric field
operators. The quantity $\bm{Q}$ within the correlation function (160) is
given by

$\bm{Q}(x,x_{1},\vec{k}_{\parallel},\omega)=\frac{\hbar\mu_{0}}{\pi}\frac{\omega^{4}}{c^{2}}\text{Im}[\varepsilon(\omega_{-})]\bm{G}(x,x_{1},\vec{k}_{\parallel},\omega)\vec{\cdot}\bm{G}^{\dagger}(x,x_{1},\vec{k}_{\parallel},\omega).$
(161)

The result (160) is that which enters Lifshitz theory (148), plus an
additional velocity dependent term equal to an integral over low frequency
excitations within the moving plate. The new term $\bm{Q}$ comes from an
imbalance between the frequency spectrum of the electrical current in the
stationary plate and the one in motion, which in turn is due to the Doppler
effect. This is the origin of the friction between the plates. The imbalance
of frequencies is similar to the effect discussed in section III.5, where the
Doppler effect causes a sign change of some of the frequencies between
reference frames, which can slow a moving particle.

Integrating the term on the second line of (160) over the moving plate and
applying the result (159) results in two terms,

$\int_{-\infty}^{0}dx_{1}\bm{Q}(x,x_{1},\vec{k}_{\parallel},\omega)=\frac{\hbar\mu_{0}\omega^{2}}{\pi}\frac{\bm{G}(x,x,\vec{k}_{\parallel},\omega)-\bm{G}^{\dagger}(x,x,\vec{k}_{\parallel},\omega)}{2i}\\\
-\frac{\hbar\mu_{0}\omega^{2}}{2\pi
i}\left\\{\bm{G}(x,0,\vec{k}_{\parallel},\omega)\vec{\cdot}\vec{e}_{x}\vec{\times}\left(\vec{\nabla}\vec{\times}\bm{G}^{\dagger}(x,0,\vec{k}_{\parallel},\omega)\right)-\text{h.c.}\right\\}$
(162)

where the curl in the second line is taken with respect to the second spatial
coordinate. The term on the second line came from an integration by parts, and
can be thought of as a flux of momentum leaving the moving plate (it takes the
form of an electric field crossed with a magnetic field, like the Poynting
vector). This is the term responsible for the frictional force between the
plates, which we now examine in detail, using the explicit form for the Green
function.

The Green function between relatively moving plates: In empty space the Green
function satisfies

$\bm{\nabla}\times\bm{\nabla}\times\bm{G}_{0}-k_{0}^{2}\bm{G}_{0}=\bm{1}\delta^{(3)}(\bm{x}-\bm{x}^{\prime})$

When written in terms of $x$, $x^{\prime}$, $\vec{k}_{\parallel}$ and $\omega$
this equals

$\bm{G}_{0}(x,x^{\prime},\vec{k}_{\parallel},\omega)=\left[\bm{1}-\frac{\vec{k}_{\pm}\otimes\vec{k}_{\pm}}{k_{0}^{2}}\right]\frac{ie^{i\sqrt{k_{0}^{2}-k_{\parallel}^{2}}|x-x^{\prime}|}}{2\sqrt{k_{0}^{2}-k_{\parallel}^{2}}}.$

where $k_{0}=\omega/c$ and

$\vec{k}_{\pm}=\text{sign}(x-x^{\prime})\vec{e}_{x}\sqrt{k_{0}^{2}-k_{\parallel}^{2}}+\vec{k}_{\parallel},$

which is simply the 1D Green function (41) with a bit of dressing to take
account of polarisation and angle of incidence. For realistic sliding
velocities (on the order of metres per second), $\bm{Q}$ is evaluated in a
regime where $k_{\parallel}\gg k_{0}$ and the field exponentially dies away
from the source333333This is known as the _electrostatic_ , or _non–retarded_
limit by those working on the physics of the electromagnetic field close to
surfaces. See, e.g. [35].. In this regime only the first few reflections (we
consider two) of the field from the two plates contribute significantly to the
Green function. Moreover when the permeability $\mu=1$ the Fresnel reflection
coefficient for $s$–polarised radiation is zero in this limit343434See e.g.
[2]. Therefore the Green function can be approximated by

$\bm{G}(x,x^{\prime},\vec{k}_{\parallel},\omega)\sim\bm{G}_{0}(x,x^{\prime},\vec{k}_{\parallel},\omega)+\frac{1}{2\kappa}\bigg{[}\vec{u}\vec{\otimes}\vec{u}^{\star}\,r_{p}(\omega_{-})e^{-\kappa(x+x^{\prime})}+\vec{u}^{\star}\vec{\otimes}\vec{u}\,r_{p}(\omega)e^{-\kappa(2a-x-x^{\prime})}\\\\[5.0pt]
+\vec{u}^{\star}\vec{\otimes}\vec{u}^{\star}\,r_{p}(\omega)r_{p}(\omega_{-})e^{-\kappa(2a+x^{\prime}-x)}+\vec{u}\vec{\otimes}\vec{u}\,r_{p}(\omega)r_{p}(\omega_{-})e^{-\kappa(2a+x-x^{\prime})}\bigg{]},\\\
$ (163)

where $\kappa=\sqrt{k_{\parallel}^{2}-k_{0}^{2}}$, $r_{p}(\omega)$ is the
$p$–polarised reflection coefficient for a body at rest, and the $p$–polarised
and $s$–polarised unit vectors are respectively given by

$\displaystyle\vec{u}$
$\displaystyle=\frac{1}{k_{0}k_{\parallel}}[k_{\parallel}^{2}\vec{e}_{x}-i\kappa\vec{k}_{\parallel}]$
$\displaystyle\vec{v}$
$\displaystyle=\frac{1}{k_{\parallel}}\vec{e}_{x}\vec{\times}\vec{k}_{\parallel}.$
(164)

The curl of the Hermitian conjugate of the Green function (LABEL:sarh:gfea)
with respect to the second index is

$\vec{\nabla}^{\prime}\vec{\times}\bm{G}^{\dagger}(x,x^{\prime},\vec{k}_{\parallel},\omega)\sim\vec{\nabla}^{\prime}\vec{\times}\bm{G}_{0}(x,x^{\prime},\vec{k}_{\parallel},\omega)-\frac{ik_{0}}{2\kappa}\vec{v}\otimes\bigg{[}\vec{u}^{\star}r_{p}^{\star}(\omega_{-})e^{-\kappa(x+x^{\prime})}+\vec{u}r_{p}^{\star}(\omega)e^{-\kappa(2a-x-x^{\prime})}\\\\[5.0pt]
+\vec{u}r_{p}^{\star}(\omega)r_{p}^{\star}(\omega_{-})e^{-\kappa(2a+x^{\prime}-x)}+\vec{u}^{\star}r_{p}^{\star}(\omega)r_{p}^{\star}(\omega_{-})e^{-\kappa(2a+x-x^{\prime})}\bigg{]}.$
(165)

Inserting expressions (LABEL:sarh:gfea) and (165) for the Green function into
the surface term in (162) we find

$\frac{1}{2i}\vec{e}_{x}\vec{\cdot}\left[\bm{G}(x,0,\vec{k}_{\parallel},\omega)\vec{\cdot}\vec{e}_{x}\vec{\times}(\vec{\nabla}\vec{\times}\bm{G}^{\dagger}(x,0,\vec{k}_{\parallel},\omega))-\text{h.c.}\right]\vec{\cdot}\vec{e}_{y}\\\
=\frac{k_{y}e^{-2\kappa
a}}{k_{0}^{2}}\text{Im}[r_{p}(\omega)]\text{Im}[r_{p}(\omega_{-})]$ (166)

where terms of order $\exp(-4\kappa a)$ and smaller have been neglected and
only the symmetric part of the tensor has been retained. The $x$–$y$ component
of the electric field correlation function (160) is

$\langle
0|\hat{E}_{x}(\vec{x},t)\hat{E}_{y}(\vec{x},t)|0\rangle=\frac{\hbar\mu_{0}}{\pi}\int\frac{d^{2}\vec{k}_{\parallel}}{(2\pi)^{2}}\int_{Vk_{y}}^{\infty}d\omega\omega^{2}\text{Im}[\bar{G}^{(S)}_{xy}(x,x,\vec{k}_{\parallel},\omega)]+\\\
\frac{\hbar\mu_{0}}{2\pi
i}\int\frac{d^{2}\vec{k}_{\parallel}}{(2\pi)^{2}}\int_{0}^{Vk{y}}d\omega\omega^{2}\vec{e}_{x}\vec{\cdot}\left[\bm{G}(x,0,\vec{k}_{\parallel},\omega)\vec{\cdot}\vec{e}_{x}\vec{\times}\left(\vec{\nabla}\vec{\times}\bm{G}^{\dagger}(x,0,\vec{k}_{\parallel},\omega)\right)-\text{h.c.}\right]\vec{\cdot}\vec{e}_{y},$
(167)

where only the symmetric part of the tensor on the second line is included. In
this ‘electrostatic’ regime the magnetic correlation function does not
contribute to the stress353535See [32].. Furthermore it was shown in [37] that
the first term on the right of (167) is zero, although we shall not prove this
here. Therefore the frictional force is $\varepsilon_{0}$ times the second
term in (167) which after we apply (166) is

$\langle
0|\hat{T}_{xy}|0\rangle=\frac{\hbar}{\pi}\int\frac{d^{2}\vec{k}_{\parallel}}{(2\pi)^{2}}\int_{0}^{Vk_{y}}d\omega
k_{y}\text{Im}[r_{p}(\omega)]\text{Im}[r_{p}(\omega_{-})]e^{-2k_{\parallel}a},$
(168)

where we have taken the limit $k_{\parallel}\gg k_{0}$. The above off diagonal
element of the stress tensor is the limiting expression found by Volokitin and
Persson [33] and Pendry [36]. As we already mentioned, this limit is
appropriate for low sliding velocities (relative to the speed of light), and
we can see that the frictional force is determined by the imaginary parts of
the reflection coefficients at low frequencies and large wave–vectors. The
integral over $\omega$ is such that $\omega_{-}\leq 0$ and therefore
$\text{Im}[r_{p}(\omega_{-})]\leq 0$ and the stress is negative. This means
that the force on the left hand body acts against its motion, and the force on
the right hand body acts in the opposite direction; i.e. the relative motion
of the bodies is being reduced. Roughly speaking, for materials with a large
degree of dissipation this frictional force is large. The vacuum field close
to a dissipative body can thus be thought of as being like a viscous fluid
that inhibits motion parallel to the surface.

##### Exercise:

Show that the reflection coefficient for $p$–polarised radiation

$r_{p}=\frac{\varepsilon(\omega)\sqrt{k_{0}^{2}-k_{\parallel}^{2}}-\sqrt{\varepsilon(\omega)k_{0}^{2}-k_{\parallel}^{2}}}{\varepsilon(\omega)\sqrt{k_{0}^{2}+k_{\parallel}^{2}}+\sqrt{\varepsilon(\omega)k_{0}^{2}-k_{\parallel}^{2}}}$

approaches $(\varepsilon(\omega)-1)/(\varepsilon(\omega)+1)$ as
$k_{\parallel}/k_{0}\to\infty$. Using this result show that in the same limit

$\text{Im}[r_{p}(\omega)]\to\frac{2\text{Im}[\varepsilon(\omega)]}{|\varepsilon(\omega)+1|^{2}}.$

For the simple case of a constant conductivity $\sigma$ the permittivity is
$\varepsilon(\omega)=\varepsilon_{b}+i\sigma/\omega\varepsilon_{0}$
($\varepsilon_{b}$ assumed constant). For this case show that the stress
tensor (168) is given by

$\langle
0|\hat{T}_{xy}|0\rangle=\frac{4\sigma^{2}\hbar}{\pi\varepsilon_{0}^{2}}\int\frac{d^{2}\vec{k}_{\parallel}}{(2\pi)^{2}}\int_{0}^{Vk_{y}}d\omega
k_{y}\frac{\omega\omega_{-}}{|i\sigma/\varepsilon_{0}+(1+\varepsilon_{b})\omega|^{2}|i\sigma/\varepsilon_{0}+(1+\varepsilon_{b})\omega_{-}|^{2}}e^{-2k_{\parallel}a}.$

Using your favourite software or programming language numerically evaluate
this integral as a function of $V$ and $a$ and plot the results.

We complete this tutorial having developed the theory of quantum friction and
quantum forces, all selfconsistently within the formalism of macroscopic QED.
This formalism is a very general way of treating the quantum mechanics of
macroscopic bodies and the electromagnetic field and may be applied to any
problem in Casimir physics. Not only does this justify the formulae of
Lifshitz theory on the basis of a complete quantum mechanical theory, but it
opens up the possibility of exploring new effects that might arise, for
example, when the centre of mass of a macroscopic body is prepared in a
quantum mechanical state.

##  IV Problems

##### Problem:

Taking the Hamiltonian (141) show that in the limit of a very massive body in
uniform motion it becomes

$H=-\vec{V}\vec{\cdot}\int
d^{3}\vec{x}\int_{0}^{\infty}d\omega\left[(\vec{\nabla}\vec{\otimes}\vec{X}_{\omega})\vec{\cdot}\vec{\Pi}_{\vec{X}_{\omega}}-\alpha(\omega,\vec{x}-\vec{R})\vec{X}_{\omega}\vec{\times}\vec{B}\right]+H_{0}+\text{const.}$
(169)

where $\vec{V}=\vec{p}/M$. Show that if the constant is neglected then this
Hamiltonian can have an arbitrarily low energy (which could become negative),
while the original Hamiltonian (141) is always positive. What is the
resolution of this apparent contradicition?

##### Problem:

Find the equations of motion for the electromagnetic field and the reservoir
from the Hamiltonian (169), and find expressions for the electromagnetic field
in terms of the $\vec{C}_{\omega}(k_{x},y,z)$ and
$\vec{C}_{\omega}^{\star}(k_{x},y,z)$ of the reservoir. Assume that the medium
is homogeneous and moves along the $x$–axis.

##### Problem:

Extend the solutions of the above problem to the quantum case and show that a
suitable Hamiltonian to describe the evolution of the operators is

$\hat{H}=\int\frac{dk_{x}}{2\pi}\int_{0}^{\infty}d\omega\hbar(\omega+Vk_{x})\hat{\vec{C}}_{\omega}^{\dagger}(k_{x},y,z)\vec{\cdot}\hat{\vec{C}}_{\omega}(k_{x},y,z).$

What are the eigenvalues of this Hamiltonian and how do they differ from the
$V=0$ case?

##### Problem:

Consider a detector coupled to the electromagnetic field with a Hamiltonian

$\hat{H}_{D}=\hbar\omega_{0}\left[\hat{a}^{\dagger}\hat{a}+\frac{1}{2}\right]$

and an interaction

$\hat{H}_{I}=i\beta\left[\hat{a}-\hat{a}^{\dagger}\right]\vec{e}\vec{\cdot}\hat{\vec{E}}(\vec{x}_{0},t),$

where $\beta$ is a constant determining the strength of the interaction, and
$\vec{e}$ is a unit vector determining its orientation. Suppose that the
detector is initially in the first excited state and find an expression (valid
to first order in perturbation theory) for the transition rate into the ground
state in terms of the electric field operator. Using the expression for the
electric field operator given in the text evaluate this in terms of the Green
function.

##### Problem:

Find an expression for the electromagnetic Green function in a region of space
outside of a dielectric half space ($x>0$) characterised by reflection
coefficients $r_{\lambda}(\omega,\vec{k})$ for the two polarisations
$\lambda=1,2$. Use the result of problem (IV) to find the transition rate of
the detector as a function of the distance from the surface.
_Hint_ : Start from the free space Green function, decompose it into a sum of
freely propagating waves, and add in the reflected waves necessary to fulfil
the boundary conditions.

##### Problem:

Find an expression for the electromagnetic Green function in a region of space
outside of a dielectric cylinder of radius $R$ characterised by some
$\varepsilon(\omega)$. Use the result of problem (IV) to find the transition
rate of the detector as a function of the distance from the cylinder.

## References

  * [1] E. M. Lifshitz and L. P. Pitaevskii, Statistical Physics (Part 2), Vol. IX (Butterworth Heinemann, Oxford, 2003).
  * [2] L. D. Landau, E. M. Lifshitz, and L. P. Pitaevskii, Electrodynamics of Continuous Media, Vol. VIII (Butterworth Heinemann, Oxford, 2004).
  * [3] J. J. Hopfield, “Theory of the contribution of Excitons to the Complex Dielectric Constant of Crystals” Phys. Rev., 112, 1555 (1958).
  * [4] B. Huttner and S. M. Barnett, “Quantization of the electromagnetic field in dielectrics” Phys. Rev. A, 46, 4306 (1992).
  * [5] L. G. Suttorp and A. J. van Wonderen, “Fano diagonalization of a polariton model for an inhomogeneous absorptive dielectric” Euro. Phys. Lett., 67, 766 (2004).
  * [6] F. Kheirandish and M. Soltani, “Extension of the Huttner–Barnett model to a magnetodielectric medium” Phys. Rev. A, 78, 012102 (2008).
  * [7] S. Scheel and S. Y. Buhmann, “Macroscopic quantum electrodynamics - Concepts and applications” Act. Phys. Solv., 58, 675 (2008).
  * [8] T. G. Philbin, “Canonical quantization of macroscopic electromagnetism” New J. Phys., 12, 123008 (2010).
  * [9] L. D. Landau and E.M. Lifshitz, Mechanics, Vol. I (Butterworth Heinemann, Oxford, 2004).
  * [10] C. Lanczos, The Variational Principles of Mechanics (Dover Publications, New York, 2013).
  * [11] L. D. Landau and E. M. Lifshitz, Statistical Physics (Part 1), Vol. V (Butterworth Heinemann, Oxford, 2005).
  * [12] L. D. Landau and E.M. Lifshitz, The Classical Theory of Fields, Vol. II (Butterworth Heinemann, Oxford, 2003).
  * [13] L. H. Ryder, Quantum Field Theory (Cambridge University Press, Cambridge, 2003).
  * [14] H. Kleinert, Path Integrals in Quantum Mechanics, Statistics, Polymer Physics and Financial Markets (World Scientific, Singapore, 2006).
  * [15] P. A. M. Dirac, Quantum Mechanics (Oxford University Press, London, 1958).
  * [16] W. Eckhardt, “First and second fluctuation–dissipation–theorem in electromagnetic fluctuation theory” Opt. Comm., 41, 305 (1982).
  * [17] U. Fano, “Effects of Configuration Interaction on Intensities and Phase Shifts” Physical Review, 124, 1866 (1961).
  * [18] S. A. R. Horsley, “Canonical quantization of the electromagnetic field interacting with a moving dielectric medium” Phys. Rev. A, 86, 023830 (2012).
  * [19] J. W. Dettman, Applied Complex Variables (Dover Publications, New York, 1984).
  * [20] S. A. R. Horsley and T. G. Philbin, “Canonical quantization of electromagnetism in spatially dispersive media” New J. Phys., 16, 013030 (2014).
  * [21] S. Weinberg, The Quantum Theory of Fields (Vol. 1) (Cambridge University Press, Cambridge, 1995).
  * [22] U. Leonhardt, Essential Quantum Optics, (Cambridge University Press, New York, 2010).
  * [23] S. M. Barnett, “Resolution of the Abraham–Minkowski Dilemma” Phys. Rev. Lett., 104, 070401 (2010).
  * [24] L. Novotny and B. Hecht, Principles of Nano Optics (Cambridge University Press, Cambridge, 2006).
  * [25] M. D. Semon and J. R. Taylor, “Thoughts on the magnetic vector potential” Am. J. Phys., 64, 1361 (1996).
  * [26] M. Wilkens, “Significance of Röntgen current in quantum optics: Spontaneous emission of moving atoms” Phys. Rev. A, 49, 570 (1994).
  * [27] R. Loudon, “Theory of the radiation pressure on dielectric surfaces” J. Mod. Opt., 49, 821(2002).
  * [28] C. K. Law, “Interaction between a moving mirror and radiation pressure: A Hamiltonian formulation” Phys. Rev. A, 51, 2537 (1995).
  * [29] W. M. R. Simpson, S. A. R. Horsley, and U. Leonhardt, “Divergence of Casimir stress in inhomogeneous media” Phys. Rev. A, 87, 043806 (2013).
  * [30] L. S. Levitov, “Van Der Walls’ Friction” Europhys. Lett., 8, 499 (1989).
  * [31] V. B. Braginsky and F. Y. Khalili, “Friction and fluctuations produced by the quantum ground state” Phys. Lett. A, 161, 197 (1991).
  * [32] J. B. Pendry, “Shearing the vacuum - quantum friction” J. Phys. Cond. Mat. 9, 10301 (1997).
  * [33] A. I. Volokitin and B. N. J. Persson, “Theory of the interaction forces and the radiative heat transfer between moving bodies” Phys. Rev. B, 78, 155437 (2008).
  * [34] V. B. Berestetskii, L. P. Pitaevskii, and E. M. Lifshitz, Quantum Electrodynamics, Vol. IV (Butterworth Heinemann, Oxford, 2004).
  * [35] I. D. Mayergoyz, Plasmon Resonances in Nanoparticles, (World Scientific, Singapore, 2013).
  * [36] J. B. Pendry, “Quantum friction - fact or fiction?” New J. Phys. 12, 033028 (2010).
  * [37] T. G. Philbin and U. Leonhardt, “No quantum friction between uniformly moving plates” New J. Phys. 11, 033035 (2009).


# Spatially-Resolved Temperature and Density Structures of Nearby Hii regions

Yifei Jin (金刈非) Research School for Astronomy & Astrophysics, Australian
National University, Canberra, Australia, 2611  ARC Centre of Excellence for
All Sky Astrophysics in 3 Dimensions (ASTRO 3D)  Ralph Sutherland  Research
School for Astronomy & Astrophysics, Australian National University, Canberra,
Australia, 2611  Lisa J. Kewley  Research School for Astronomy & Astrophysics,
Australian National University, Canberra, Australia, 2611  Institute for
Theory and Computation, Harvard-Smithsonian Center for Astrophysics,
Cambridge, MA 02138, USA  ARC Centre of Excellence for All Sky Astrophysics in
3 Dimensions (ASTRO 3D)  David C. Nicholls  Research School for Astronomy &
Astrophysics, Australian National University, Canberra, Australia, 2611  ARC
Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D)

(Received January 1, 2018; Revised January 7, 2018)

###### Abstract

Photoionization models frequently assume constant temperature or density
within HII regions. We investigate this assumption by measuring the detailed
temperature and density structures of four Hii regions in the Large Magellanic
Cloud and the Small Magellanic Cloud, using integral-field spectroscopic data
from the Wide-Field Spectrograph on the ANU 2.3m telescope. We analyse the
distribution of emission-lines of low-ionization species, intermediate-
ionization species and high-ionization species. We present the complex
electron temperature and density structures within Hii regions. All four
nebulae present a negative gradient in the electron density profile. Both
positive and negative temperature gradients are observed in the nebulae. We
create a series of nebula models with a constant ISM pressure and varying
temperature and density distributions. Comparison of the line ratios between
our Hii regions and models suggests that none of the simple nebula models can
reproduce the observed temperature and density structures. Comparison between
the models and the data suggests that the ISM pressure of nebulae in LMC and
SMC is between log($P/k$)=6-7.5. Complex internal structures of the nebulae
highlight the importance of future Monte-Carlo photoionization codes for
accurate nebula modeling, which include a comprehensive consideration of
arbitrary geometries of Hii regions.

ISM: HII regions — ISM: structure — galaxies:starburst

††journal: ApJ

##  1 Introduction

To fully understand galaxies, we need to accurately model the Hii regions
within galaxies. The electron temperature and density are two of the
fundamental quantities that determine the nebular emission-line spectrum
within Hii regions. Accurate predictions of the interstellar medium (ISM)
physical quantities rely on the assumed or modeled distributions of electron
temperature, $T_{e}$, and density, $n_{e}$.

The ISM electron temperature and density distributions across Hii regions are
complex. The density distributions in most Hii regions show filamentary or
shell-like structures (Kennicutt, 1984). The density measurements within these
structures have diverse radial gradients (Binette et al., 2002; Herrera-Camus
et al., 2016; Rubin et al., 2016), although relatively uniform density
distributions are seen in some cases (Garnett & Dinerstein, 2001; Krabbe &
Copetti, 2002; García-Benito et al., 2010). Surveys of both galactic and
extragalactic Hii regions indicate that approximately half the Hii regions
have internal variations in electron density (Copetti et al., 2000; Malmann et
al., 2002).

The presence of temperature fluctuations in Hii regions is under significant
debate due to the notorious difficulty in measuring the electron temperature.
Some studies show evident temperature fluctuations (Peimbert et al., 2004;
Hägele et al., 2006; Kewley et al., 2019b), while others find no temperature
fluctuations (Stasińska et al., 2013; Liu et al., 2006).

The underlying temperature and density structures in Hii regions affect the
robustness of metallicity calibrations. The most direct calibration method is
to convert metallicities from electron temperatures in the $\rm O^{2+}$ zone
after correcting for the unseen ionization stages (the so-called ionization
correction factor, ICF). The ICF is determined from the photoionization models
under the assumption of constant electron temperature and density across the
Hii region. Hägele et al. (2008) find electron temperature variations raise
the uncertainty of metallicities by up to 0.2 dex through the ICF method.

The temperature variation is one of the potential causes of the discrepancies
among metallicities measured from different emission-lines (Kewley et al.,
2019b). Commonly used metallicity-sensitive emission-line ratios are ([O
ii]$\lambda 3727$ \+ [O iii]$\lambda\lambda 4959,5007$)/H$\beta$ (known as
$R_{23}$), [O iii]$\lambda 5007$/H$\beta$/[N ii]$\lambda 6584$/H$\alpha$(known
as O3N2) and [N ii]$\lambda 6584$/H$\alpha$. Recombination lines give larger
metallicities than auroral lines. The offset between recombination line and
auroral metallicities has been suggested as arising from temperature
fluctuations in Hii regions (Peimbert, 1967), which is known as the abundance
discrepancy factor (ADF).

The ISM pressure incorporates the complex structures of electron temperature
and density. Constant pressure models are applicable to the majority of Hii
regions where the sound-crossing time is shorter than the heating and cooling
timescale. Pressure photoionization models with complex temperature and
density structures have been successfully applied to galaxies at $z\sim 2$,
suggesting the existence of the high-pressure Hii regions in high-$z$ galaxies
(Lehnert et al., 2009).

Current photoionization models significantly underestimate the complexity of
the internal structures of Hii regions. Isothermal or constant density ISM
conditions are often assumed in photoionization models. Isothermal models with
a single temperature across the nebula have been used in calibrating the
standard density-sensitive line ratios (Osterbrock, 1989; Rubin et al., 1994),
and determining the excitation power source in galaxies (Filippenko, 1985;
Groves et al., 2004). The constant density distribution assumption is adopted
in measuring the electron temperature and density in Hii regions (Osterbrock,
1989). The constant density assumption is also used in deriving the
metallicities measured from auroral lines. In modern constant pressure models,
the majority of modeled Hii regions produce complex internal gradients in
electron temperature and density, while only few specific cases with
log(O/H)+12 $\sim 8.23$ can be approximated with the isothermal or constant
density conditions (Kewley et al., 2019a).

Nearby Hii regions are laboratories for understanding the internal structure
of the ISM electron temperature and density. With the advantage of proximity,
the Magellanic Clouds provide a rich reservoir of Hii regions for better
understanding the ISM temperature and density fluctuations (Kennicutt, 1984;
Kennicutt & Hodge, 1986; Chu & Kennicutt, 1986; Vermeij & van der Hulst, 2002;
Peña-Guerrero et al., 2012b). Integral-field unit (IFU) observations offers
highly spatially-resolved data of Hii regions in the Magellanic Clouds for the
detailed study of electron temperature and density structures (Dopita &
Seitenzahl, 2019; Kumari et al., 2017). Future IFU data sets of nearby Hii
regions will be released in the on-going TYPHOON program (Seibert et al. in
prep.) and the up-coming SDSS-V/Local Volume Mapper (LVM) (Kollmeier et al.,
2017).

In this paper, we present deep IFU observations of four Hii regions in the
Large Magellanic Cloud (LMC) and the Small Magellanic Cloud (SMC). The data
were obtained using the Wide-Field Spectrograph (WiFeS) with grating
resolution of $R\sim 7000$. We derive the electron temperature and density
with the deep IFU data, and further analyze the thermal and density structures
within each Hii region using the highly spatially-resolved datacube. This is
the first study of nebular internal structures with a sample of LMC/SMC
nebulae using the IFU data with a resolution of $\rm 0.2-0.3~{}pc$. This paper
is structured as below. We first describe the sample of Hii regions in Section
2 and describe the data and the data reduction in Section 3. We describe the
calibrations of electron temperature, density and pressure in Section 4. We
present the electron temperature and density structures of Hii regions in
Section 5. In Section 6, we compare our Hii regions with the pressure
photoionization models. The discussion is given in Section 7.

##  2 The sample

Figure 1: Collapsed IFU images of four Hii regions in LMC and SMC. The field
of view is $25^{\prime\prime}\times 38^{\prime\prime}$. Grey scale color
indicates the line-of-sight integrated flux.

The four Hii regions are selected from a series of spatially-resolved
observations of the nebulae in the Magellanic Clouds (Dopita et al., 2016,
2018). The four Hii regions are bright and young compact Hii regions ionized
by a single O-star or simple star clusters. Three Hii regions, N456, N77A and
N88A are from SMC, and one, N191A, is from LMC. The simple ionizing sources
make these compact Hii regions ideal targets for testing current
photoionization models. The compact size ($\sim 10^{\prime\prime}$ in
diameter) and proximity allows these compact Hii regions to fit within the
WiFeS field of the view (FoV). Some fundamental physical quantities derived
from the global spectra of these Hii regions are listed in Table 1.

LMC N191 is a complex Hii region which presents two components in H$\alpha$,
[O iii] and H$\beta$ narrow band images (Henize, 1956). N191A is more compact
and brighter between the two components, which has a diameter of
$5^{\prime\prime}.2$ (1.2 pc at LMC distance) (Selier & Heydari-Malayeri,
2012). The broad $V$-band image reveals that the ionized gas in N191A is
excited by an O-star located at the center of the nebula (Selier & Heydari-
Malayeri, 2012).

SMC N77A is a spherical Hii region having a radius of $10^{\prime\prime}$
(Selier & Heydari-Malayeri, 2012). A dust lane crosses SMC N77A along the
east-west direction (Selier & Heydari-Malayeri, 2012; Toribio San Cipriano et
al., 2017). The major ionizing source of N77A is a B-star residing north of
the dust lane. The average extinction of N77A is $\rm A_{V}$$\sim 0.25$ mag
(Selier & Heydari-Malayeri, 2012).

SMC N88A is the brightest Hii region in the SMC (Testor et al., 2003). N88A is
a $\sim 1$-pc sized ionized gas “blob” excited by an O3-4 star (Testor et al.,
2003). The heavy dust core in N88A produces its high integral extinction of
$\rm A_{V}$$=1.38$ mag (Garnett et al., 1995; Testor et al., 2003).

SMC N456 is the second brightest Hii region in the SMC (Peimbert & Torres-
Peimbert, 1976), with an inhomogeneous thermal structure. Long-slit
spectrographic study indicates the mean-square temperature variation, $t^{2}$
(first introduced by Peimbert (1967)), in SMC N456 is $0.067\pm 0.013$ (Peña-
Guerrero et al., 2012a).

Table 1: The fundamental physics quantities of the Hii region sample aaThe data are collected from the literature Selier & Heydari-Malayeri (2012) for LMC N191A and SMC N77A, Testor et al. (2003) for SMC N88A and Peña-Guerrero et al. (2012a) for SMC N456. Hii Region | Size | $\rm A_{V}$ | $T_{e}$ | $n_{e}$
---|---|---|---|---
| pc | mag | K | cm-3
LMC N191A | 1.2 | $\sim 2.70$ | 10,800 | 600
SMC N77A | 2.9 | $\sim 0.25$ | 14,240 | 60
SMC N88A | 1.2 | $\sim 1.38$ | $\sim 14,100$ | 2870
SMC N456 | – | – | $\sim 12,000$ | $\sim 100$

##  3 Observation and data reduction

The spectroscopic data were observed using WiFeS. WiFeS (Dopita et al., 2007,
2010) is an image-slicing integral-field spectrograph mounted at the
Australian National University 2.3 m telescope at the Siding Spring
Observatory. WiFeS has a $25^{\prime\prime}\times\ 38^{\prime\prime}$ field of
the view with spaxels each $1\times\ 1^{\prime\prime}$ in angular size. The
double-beam spectrograph feeds the light into the blue and red cameras
simultaneously to cover the wavelength ranging from 3300 to 9000 Å. We used
the high-resolution gratings of $R=7000$ by splitting the spectra into $U$,
$B$, $R$ and $I$ bands.

The observations were taken between November and December 2015 with seeing
FWHM 1.0-2.0 arcsec. Each Hii region was observed using 10, 100 and 1000 (or
800) second exposures. The longer exposures were aimed at securing flux
measurements of faint emission lines down to $\sim 10^{-5}$ of the H$\beta$
flux. The shorter exposures were necessary to minimize the saturation on the
bright hydrogen and oxygen lines.

The data were reduced using the standard WiFeS pipeline, PyWiFeS (Childress et
al., 2014) . Telluric lines were removed by subtracting observed clear sky
regions adjacent to the Hii regions, from each observation. The separate
camera spectra were then combined using the overlapping spectral regions
between adjacent cameras. The final data has a spectral resolution of 43 $\rm
kms^{-1}$ and a spatial resolution of 0.3 pc for SMC N77A, SMC N88A and SMC
N456, and a spatial resolution of 0.24 pc for LMC N191A.

##  4 Derived quantities

###  4.1 Emission-line maps

Figure 2: Emission-line maps of SMC N77A. We show the maps of Balmer lines
H$\alpha$ and H$\beta$, the emission-lines of high-ionization species [O
iii]$\lambda 5007$ and [Ne iii]$\lambda 3869$, the emission-lines of
intermediate-ionization species [Ar iii]$\lambda 7136$, [S iii]$\lambda 6312$,
and the emission-lines of low-ionization species [O ii]$\lambda 3727$, [N
ii]$\lambda 6584$ and [S ii]$\lambda 6716$. Figure 3: Same as Figure 2 but for
LMC N191. Figure 4: Same as Figure 2 but for SMC N88A. Figure 5: Same as
Figure 2 but for SMC N456. The central pixel in the [O iii]$\lambda$5007 map
is a bad pixel because the [O iii]$\lambda$5007 is saturated at the center of
SMC N456.

The emission-line fluxes are obtained from each pixel by integrating the
spectrum subtracted from the underlying stellar continuum and nebular
continuum. Continuum components are estimated by using the average flux over a
10Å wavelength interval on either the left side or the right side to each
emission-line. In Figure 2 to Figure 5, we show maps of the emission-lines
representing species in the high-ionization zone, the intermediate-ionization
zone and the low-ionization zone of the nebula. We use the Balmer lines to
trace the main body of Hii regions. The [O iii]$\lambda 5007$ and [Ne
iii]$\lambda 3869$ lines are used to trace the high-ionization zone. The [Ar
iii]$\lambda 7136$ and [S iii]$\lambda 6312$ lines are the tracers of the
intermediate-ionization zone. The [O ii]$\lambda 3727$, [N ii]$\lambda 6584$
and [S ii]$\lambda 6716$ emission-lines are used to trace the low-ionization
zone of the nebula. The boundary of the Hii region is identified at the pixels
with $\rm Flux(H\alpha)=$10%$\rm Flux(H\alpha)_{peak}$, where $\rm
Flux(H\alpha)_{peak}$ is the maximal H$\alpha$ flux.

In Figure 2 to Figure 5, a clear spatial segregation exists among different
ionization zones. The high-ionization zone that is traced by [O iii]$\lambda
5007$ and [Ne iii]$\lambda 3869$ is compact around the nebular center because
of the hard ionization field from the ionizing star in proximity to the star.
In contrast, the low-ionization zone traced by [O ii]$\lambda 3727$, [N
ii]$\lambda 6584$ and [S ii]$\lambda 6716$ is located at the boundary of the
Hii regions. The [Ar iii]$\lambda 7136$ and [S iii]$\lambda 6312$ lines are
located between the high-ionization and low-ionization zones.

Our Hii region sample presents diverse nebular morphologies. The emission-
lines in LMC N191 have an extended morphology, which is consistent with the
early H$\alpha$ image observed by Henize (1956). SMC N88A appears as a bright
compact point source for the emission-lines of high-ionization and
intermediate-ionization species, which is consistent with the near-infrared
imaging of SMC N88A (Testor et al., 2010). For the emission-lines of low-
ionization species, an extended region exists north-east of the bright nebular
center. SMC N456 presents a shell-like morphology for the emission-lines of
low-ionization species. However, the emission-lines of high-ionization species
present a spherical morphology. The extended morphology of SMC N456 is also
confirmed in Arellano-Córdova & Rodríguez (2020). SMC N77A appears to have a
two-lobe morphology for emission-lines of high-ionization and intermediate-
ionization species, which is split by a dust lane crossing from east to west
(Selier & Heydari-Malayeri, 2012). For the low-ionization species, the
emission-lines are located on the nebular edge, leading to a shell-like
morphology.

###  4.2 Extinction Correction

Figure 6: Maps of attenuation $\rm A_{V}$. The color code indicates the value
of $\rm A_{V}$. The white contours present the distribution of H$\beta$
emission-line fluxes.

The extinction correction is carried out through the H$\alpha$/H$\beta$ Balmer
decrement, using the Fitzpatrick (1999) extinction law. We select a Rv$=3.16$
for nebula in LMC and Rv$=2.93$ for nebula in SMC (Pei, 1992). The Fitzpatrick
extinction curve is a generally-used parameterized curve. Compared to the
average extinction curves of the LMC and the SMC (Gordon et al., 2003), we
find the difference of the extinction corrected fluxes is 0.7%-8% by applying
different extinction curves on our sample.

Figure 6 presents maps of the derived $\rm A_{V}$ with the distribution of
H$\beta$ fluxes. We notice that the enhancement of $\rm A_{V}$$\sim 0.7$ $mag$
in the SMC N77A is at the position of the central dust lane. In SMC N88A, the
maximum $\rm A_{V}$ is $\sim$2.5 $mag$, located at the peak of H$\beta$
fluxes.

###  4.3 Determining electron density

We measure the electron density, $n_{e}$, for each resolution element by using
the calibrations in Osterbrock & Ferland (2006) (hereafter OF06) based on
single-atom fixed temperature, fixed density models, and the diagnostics given
by Kewley et al. (2019a) (hereafter K19) based on their theoretical density
models. Compared with single-atom models, the density models in Kewley et al.
(2019a) are more representative to the physical conditions in nebulae in terms
of the temperature structures changing with the ionization field and the
metallicity. The electron densities derived from the Osterbrock & Ferland
(2006) method and the Kewley et al. (2019a) density models are independent
from the ISM pressure, which can be used to compare with the constant pressure
models where the ISM pressure is a free parameter.

The electron densities are derived from the [S ii]$\lambda 6716$/ [S
ii]$\lambda 6731$ ratio and the [O ii]$\lambda 3726$/ [O ii]$\lambda 3729$
ratio respectively. We use the amplitude-to-noise (A/N) ratio to describe the
quality of each emission-line. The A/N ratio is the ratio of the peak of
emission-line flux to the noise of emission-lines, tracing both the signal-to-
noise ratio and the strength of emission-line fluxes. The A/N ratio is
sensitive to the emission-lines with low signal-to-noise ratios but evident
excess of fluxes superimposed on the continuum. We require A/N of the [O ii]
doublets and the [S ii] doublets larger than 5 in each spaxel.

The OF06 models use the atomic data from Mendoza (1983), which provides a
compilation of transition probabilities, electron excitation coefficients and
photoionization cross sections from sophisticated computation and
measurements. The K19 models use the CHIANTI v.8.0 atomic data (Del Zanna et
al., 2015) in conjunction with the energy level information of each ions from
NIST database. The K19 models are suitable for calculating the electron
density between 1 and $\rm 10^{8}~{}cm^{-3}$.

###  4.4 Determining electron temperature

We measure the electron temperatures, $T_{e}$, from collisional lines using
the method described in Osterbrock & Ferland (2006) and the calibration in
Nicholls et al. (2020) (hereafter NKS20). The Osterbrock & Ferland (2006)
method is a simple polynomial fit to the single-atom models within a
temperature range from 5000 K to 25,000 K. Nicholls et al. (2020) re-calibrate
the relationship between the [O iii] emission-line ratios and $T_{e}$ by using
a rational polynomial as shown in equation 1

${\rm
log}_{10}(T_{e})=\frac{P_{0}+P_{1}x+P_{2}x^{2}}{1+Q_{1}x+Q_{2}x^{2}+Q_{3}x^{3}},$
(1)

where $x$ is $\rm log_{10}$([O iii]$\lambda 4363$/[O iii]$\lambda\lambda
4959,5007$), and $P_{0}$, $P_{1}$, $P_{2}$, $Q_{1}$, $Q_{2}$, $Q_{3}$ are the
fit coefficients. We require A/N$>$5 for [O iii]$\lambda\lambda$4959,5007 and
[O iii]$\lambda$4363

The OF06 models use the atomic data from Mendoza (1983). The NKS20 models use
the CHIANTI v.8.0 atomic data (Del Zanna et al., 2015) in conjunction with the
energy level information of each ions from NIST database.

###  4.5 Determining $P/k$

The pressure within Hii regions is derived from the temperature, $T$, and the
particle density, $n$, of nebula through the ideal gas assumption $P=nkT$. We
assume that $T\equiv T_{e}$ and approximate the total particle density $n$
from the electron density $n_{e}$ through $n\sim 2n_{e}(1+He/H)$ given
hydrogen and helium are the majority species in Hii regions. Finally, the
nebular pressure is derived through $P/k=2n_{e}(1+He/H)T_{e}$.

##  5 Analysis

Figure 7: Maps of electron density $n_{e}$. The color code indicates the value
of $n_{e}$. First row: Maps of the electron density derived from the [S ii]
ratio based on the method in Osterbrock & Ferland (2006). Second row: Maps of
errors of the electron density derived from the [S ii] ratio based on the
method in Osterbrock & Ferland (2006). Third row: Maps of the electron density
derived from the [S ii] ratio based on the method in Kewley et al. (2019a).
Fourth row: Maps of errors of the electron density derived from the [S ii]
ratio based on the method in Kewley et al. (2019a). Figure 8: Maps of electron
density $n_{e}$. The color code indicates the value of $n_{e}$. First row:
Maps of the electron density derived from the [O ii] ratio based on the method
in Osterbrock & Ferland (2006). Second row: Maps of errors of the electron
density derived from the [O ii] ratio based on the method in Osterbrock &
Ferland (2006). Third row: Maps of the electron density derived from the [O
ii] ratio based on the method in Kewley et al. (2019a). Fourth row: Maps of
errors of the electron density derived from the [O ii] ratio based on the
method in Kewley et al. (2019a). Figure 9: Maps of electron temperature
$T_{e}$. The color code indicates the value of $T_{e}$. First row: Maps of
electron temperature derived from the [O iii] ratio based on the method in
Osterbrock & Ferland (2006). Second row: Maps of the error of electron
temperature measured from the [O iii] ratio based on the method in Osterbrock
& Ferland (2006). Third row: Maps of electron temperature derived from the [O
iii] ratio based on the method in Nicholls et al. (2020). Fourth row: Maps of
the error of electron temperature measured from the [O iii] ratio based on the
method in Nicholls et al. (2020).

###  5.1 Electron density structures

In Figure 7 and Figure 8, we present the maps of electron density derived from
the [O ii] ratio, $n_{e}(\rm OII)$, and the [S ii] ratio, $n_{e}(\rm SII)$. In
LMC N191, SMC N77A and SMC N88A, the 2D density distribution has a dense core
at the center of the nebula and the density decreases along the radius of
nebula. In SMC N456, there are three dense clumps located at three corners of
the nebula.

We measure the gradients of the electron density profiles using the following
function,

$\frac{n_{e}}{cm^{-3}}=\alpha+\beta\left(\frac{r}{pc}\right),$ (2)

where $\beta$ is the gradient of the profile. The center of nebula is
identified as the position of the peak flux of H$\alpha$. The best-fit
parameters are listed in Table 3.

The derived density gradient changes by using different methods of the
measurement. LMC N191 and SMC N88A have evident negative density gradients
with a slope of -75 $\rm cm^{-3}~{}pc^{-1}$$<\beta<$-45$\rm cm^{-3}~{}pc^{-1}$
for LMC N191 and -850 $\rm cm^{-3}~{}pc^{-1}$$<\beta<$-175 $\rm
cm^{-3}~{}pc^{-1}$ for SMC N88A. In LMC N191, the density gradient derived
from the K19 [O ii]-ratio method is smaller than the gradient derived from
other methods. In contrast, the SMC N77A and SMC N456 have the shallow and
even flat density gradients. In SMC N77A, the density gradient derived from
the [S ii]-ratio is marginally negative by -8$\pm$5 $\rm
cm^{-3}~{}pc^{-1}$$<\beta<$-6$\pm$5 $\rm cm^{-3}~{}pc^{-1}$. The gradient
derived from the [O ii]-ratio is flat with a $\beta\sim$-4$\pm$5.5 $\rm
cm^{-3}~{}pc^{-1}$. In SMC N456, the density gradient is flat within the error
range by $\beta\sim$-3$\pm$3 $\rm cm^{-3}~{}pc^{-1}$.

Table 2: The average and standard deviation of the temperature by pixels of
each Hii region
(K) | LMCN191 | SMCN77A | SMCN88A | SMCN456
---|---|---|---|---
$<T_{e}>_{OF06}$ | 8864 | 11405 | 12932 | 11321
$\sigma(T_{e})_{OF06}$ | 813 | 630 | 1763 | 423
$<T_{e}>_{NKS20}>$ | 8923 | 11490 | 13030 | 11405
$\sigma(T_{e})_{NKS20}$ | 823 | 636 | 1779 | 427

###  5.2 Temperature structures

In Figure 9, we present the maps of electron temperatures derived from [O
iii]-ratio through the methods given by Osterbrock & Ferland (2006) and
(Nicholls et al., 2020). The temperature is evenly distributed across the
entire nebula for the four Hii regions.

We calculate the average and the standard deviation of $T_{e}$ of each Hii
region. As shown in Table 2, the average temperature of the Hii region is
around 8 000 K to 13 000 K, which is consistent with the previous measurement
from the integrated spectra (Testor et al., 2003; Peña-Guerrero et al., 2012a;
Selier & Heydari-Malayeri, 2012). In each nebula, the measured temperature is
consistent between the OF06 method and the NKS20 method.

We derive the temperature gradient of each Hii region by fitting temperature
radial profiles with an analytical function as below,

$\frac{T_{e}}{K}=\alpha+\beta\left(\frac{r}{pc}\right),$ (3)

where $r$ is the radius of nebula. The center of nebula is identified as the
position of the peak flux of H$\alpha$. The best-fit parameters are listed in
table 3.

The temperature gradient is flat in the four Hii regions. The gradient is
negligible to the average temperature of each nebula by $\beta$ $\approx$
5%$<T_{e}>$, where $<T_{e}>$ is the measured average temperature. LMC N191,
SMC N77A and SMC N88A show a positive trend in the temperature radial profile
with the gradient ranging from 37$\pm$41 K $\rm pc^{-1}$ to 180$\pm$320 K $\rm
pc^{-1}$. SMC N456 shows a negative trend in the temperature radial profile
with a gradient of -27$\pm$30 K $\rm pc^{-1}$. The gradient values are smaller
than the fitting errors in the four nebulae.

Table 3: Parameters of the best-fits to the radial profiles of electron
temperature and electron densitya.

| LMC N191 |  | SMC N77A |  | SMC N88A |  | SMC N456
---|---|---|---|---|---|---|---
| $\alpha$ | $\beta$ |  | $\alpha$ | $\beta$ |  | $\alpha$ | $\beta$ |  | $\alpha$ | $\beta$
$T_{e}(OIII)$ |  |  |  |  |  |  |  |  |  |  |
(OF06) | 8570$\pm$256 | 73$\pm$87 |  | 11117$\pm$156 | 37$\pm$41 |  | 12299$\pm$7 | 179$\pm$317 |  | 11610$\pm$252 | -27$\pm$30
$T_{e}(OIII)$ |  |  |  |  |  |  |  |  |  |  |
(NKS20) | 8626$\pm$259 | 75$\pm$88 |  | 11200$\pm$157 | 37$\pm$41 |  | 12391$\pm$8 | 181$\pm$320 |  | 11697$\pm$253 | -27$\pm$31
$n_{e}(SII)$ |  |  |  |  |  |  |  |  |  |  |
(OF06) | 582$\pm$21 | -54$\pm$5 |  | 143$\pm$49 | -8$\pm$5 |  | 2375$\pm$179 | -396$\pm$117 |  | 156$\pm$45 | -3$\pm$3
$n_{e}(SII)$ |  |  |  |  |  |  |  |  |  |  |
(K19) | 748$\pm$17 | -72$\pm$6 |  | 88$\pm$51 | -6$\pm$5 |  | 1276$\pm$96 | -172$\pm$40 |  | 112$\pm$52 | -3$\pm$3
$n_{e}(OII)$ |  |  |  |  |  |  |  |  |  |  |
(OF06) | 579$\pm$13 | -46$\pm$5 |  | 107$\pm$37 | -4$\pm$5 |  | 3033$\pm$717 | -722$\pm$353 |  | 115$\pm$30 | -3$\pm$3
$n_{e}(OII)$ |  |  |  |  |  |  |  |  |  |  |
(K19) | 655$\pm$2 | -5$\pm$2 |  | 73$\pm$54 | -4$\pm$6 |  | 3244$\pm$483 | -831$\pm$189 |  | 61$\pm$38 | -2$\pm$4

a The radial profiles are fitted with the analytical function
$\rm\frac{y}{K/cm^{-3}}=\alpha+\beta\left(\frac{r}{pc}\right)$, where $y$ is
the electron temperature or the electron density, $\alpha$ is the intercept of
the radial profile and $\beta$ is the gradient of the radial profile.

##  6 Comparison to models

We combined stellar atmosphere models with the MAPPINGS V photoionization code
to create a set of Hii region models. The modeled Hii regions have a broad
range of ionization parameter and pressure.

###  6.1 Stellar ionizing source

Ionizing sources of each Hii region are identified through the photometry of
stars. In LMC N191, the dominant star is a $M_{V}=-5.27$ $mag$ O-star, of
which the spectral type is classified as O5 V (Selier & Heydari-Malayeri,
2012). The ionizing source in SMC N77A is suggested as an early-type B-star
with an intrinsic color of B-V=-0.28 (Selier & Heydari-Malayeri, 2012). Testor
et al. (2003) suggest that the ionizing source in SMC N88A has a spectral type
of O6-O8 V star with a $M_{V}$ of -4.85. The stellar context in SMC N456 is
complex because it contains four bright O-type stars, with $M_{V}$ ranging
from -4.5 to -4.75 $mag$ (Testor & Lortet, 1987).

We use the ionizing spectra from the TLUSTY stellar atmosphere library (Hubeny
& Lanz, 1995, 2017) for our photoionization models. The TLUSTY models are
fully consistent, non-LTE (local thermodynamic equilibrium) line-blanketed
atmosphere models allowing a broad range of effective temperature from 10 K to
$10^{8}$ K. We use the stellar atmosphere models with [Fe/H] abundance of -1.0
for the SMC nebulae and -0.7 for the LMC nebulae (Cioni, 2009).

###  6.2 Photoionization model grid

We used the MAPPINGS V photoionization code (Binette, 1985; Sutherland &
Dopita, 1993; Dopita et al., 2013) for modeling the ISM in compact Hii
regions. MAPPINGS V code utilizes the latest CHIANTI v.8.0 atomic database
(Del Zanna et al., 2015) for the lightest 30 elements.

MAPPINGS V is a photoionization code which is used to calculate the thermal
and ionization structures of the ISM plasma with the equilibrium of cooling
and heating, excitation and de-excitation. MAPPINGS V generates the nebular
radial variation of the electron density, temperature, the ionic fraction and
the emissions of 30 elements, under the assumption of the ISM density and
pressure structure across the nebula. MAPPINGS calculates the radiation from
two components, the ionizing source and the ionizing recombination emission.
The fractional contribution of each component is determined based on the
nebular geometry.

We construct the photoionization model grid with a set of constant ISM
pressures ranging from $\log(P/k)$=5 to $\log(P/K)$=9 in a step of 0.05 dex, a
set of ionization parameters from $\log(Q)$=6.0 to $\log(Q)$=9.5 in a step of
0.125 dex, and a varying stellar effective temperatures from $T_{\rm eff}$=30
000 K to $T_{\rm eff}$=50 000 K in a step of 2500 K. Compared with the
isochoric model which assumes a constant ISM density, constant ISM pressure is
a more realistic assumption to the majority of Hii regions where the heating
and cooling timescales are longer than the sound-crossing time (Field, 1965;
Begelman & McKee, 1990).

We use the spherical geometry models for our comparison. MAPPINGS is designed
to generate the plane-parallel and spherical nebular models. The plane-
parallel models are usually applied to the nebula whose radius is far larger
than its thickness, while the spherical models are used for the cloud powered
by a point source at its center. Our observed four Hii regions have a round
shape, with both LMC N191 and SMC N88A appearing spherical. Therefore the
spherical geometry is more appropriate than plane parallel geometry for our
HII regions.

The ISM in the Magellanic Clouds have a non-solar element ratio which cannot
be simply derived by scaling the solar metallicity. We used the measured ISM
metallicity in Magellanic Clouds from Russell & Dopita (1992). In general, the
mean metallicity of the LMC is 0.2 dex lower than our Galactic ISM, while the
SMC has a metallicity 0.6 dex lower. We included 15 elements in our LMC model
and 16 elements in our SMC model. The details of element abundances are listed
in Table 4

Table 4: The abundance of the LMC and the SMC from Russell & Dopita (1992) Element | LMC | SMC
---|---|---
| log(X/H) | log(X/H)
H | 0.0 | 0.0
He | -1.05 | -1.09
C | -3.96 | -4.24
N | -4.86 | -5.37
O | -3.65 | -3.97
Ne | -4.39 | -4.73
Na | -4.85 | -5.92
Mg | -4.53 | -5.01
Al | – | -5.60
Si | -4.28 | -4.69
S | -5.29 | -5.41
Cl | -7.23 | -7.30
Ar | -5.71 | -6.29
Ca | -6.03 | -6.16
Fe | -4.77 | -5.11
Ni | -6.04 | -6.14

###  6.3 Iterative Searching process

We take a three-step fitting process to determine the best-fit model of the
global emission-line spectrum of each Hii region. The global spectrum is
integrated across the entire field-of-view. First, we use a coarse model grid
to fit the temperature-sensitive line ratios, [N ii]$\lambda\lambda$6548,84/[N
ii]$\lambda$5755, [O iii]$\lambda$5007/[O iii]$\lambda$4363, to roughly
estimate the temperature of nebula and derive the theoretical full hydrogen
spectrum based on the temperature. The selected temperature-sensitive line
ratios are free from dust extinction. The [N ii]$\lambda\lambda$6548,84/[N
ii]$\lambda$5755 ratio traces the temperature of the low-ionization zone of
the nebula. The [O iii]$\lambda$5007/[O iii]$\lambda$4363 ratio traces the
temperature of the high-ionization zone of the nebula. The coarse model grid
has 5$<logP/K<$ 9 in a step of 0.5dex, 6$<logQ<$9.5 in 0.5dex steps and 30 000
K$<T_{\rm eff}<$50 000 K in a step of 25000 K.

Then, we perform the extinction correction by comparing the full hydrogen
spectrum between the selected model and the observation. The extinction
correction based on the Balmer decrement of H$\alpha$/H$\beta$=2.86, needs to
assume the nebula is optically thick with a constant temperature of $10^{4}$ K
(the Case-B nebula). Instead, we fit the full hydrogen spectrum by using the
Hydrogen emission-lines simultaneously, which is free from the arbitrary Case
B assumption. The observed fluxes and the best-fit model fluxes of the
hydrogen spectrum are listed in Table 5 to Table 8.

The extinction correction is taken by converting the intensity of emission-
lines, $I(\lambda)$, to the intrinsic intensity, $I_{0}(\lambda)$, through the
equation below.

$I_{0}(\lambda)=I(\lambda)\times 10^{C_{n}\times f(\lambda)},$ (4)

where $C_{n}$ is the logarithmic nebula reddening constant and $A(\lambda)/AV$
is the attenuation curve proposed by (Fitzpatrick, 1999). $f(\lambda)$ is the
extinction function normalized by H$\beta$ $f(\lambda)$ =
$A(\lambda)/AV)/(Av/A(H_{\beta})-1$, to avoid the difference between the
stellar extinction function and the nebular extinction function.

Finally, we search for the best-fit model of the de-reddened emission-line
spectrum by using the fine model grid described in Section 6.2. The best-fit
model is identified with the minimal $\chi^{2}$ value, where $\chi^{2}$ is
defined as below,

$\chi^{2}=\sqrt{\sum_{\lambda}\left(\frac{F_{\lambda}^{Obs}-F_{\lambda}^{Mod}}{F_{\lambda}^{Obs}}\right)^{2}},$
(5)

where $F_{\lambda}^{Obs}$ and $F_{\lambda}^{Mod}$ are observed and modeled
fluxes used at each step.

Table 5: LMC N191 Hydrogen Spectrum De-Reddened and Residuals Wavelength | Obs. Flux | Error | De–Red | De–Red |  $|$ Diff. $|$
---|---|---|---|---|---
(Å) | (H$\beta$= 100) | 1$\sigma$ % | Obs. Flux | Model | %
3734.375 | 2.427 | 0.91 | 2.740 | 2.397 | 12.52
3750.158 | 3.046 | 1.39 | 3.433 | 3.048 | 11.21
3770.637 | 3.775 | 1.43 | 4.246 | 3.963 | 6.66
3835.391 | 6.563 | 1.56 | 7.330 | 7.292 | 0.52
3970.079 | 14.89 | 0.79 | 16.40 | 15.86 | 3.29
4101.742 | 24.42 | 0.65 | 26.53 | 25.83 | 2.64
4340.471 | 44.03 | 0.19 | 46.68 | 46.74 | 0.13
4861.333 | 100.00 | 0.94 | 100.00 | 100.00 | 0.00
6562.819 | 331.05 | 0.23 | 285.02 | 287.73 | 0.95
8437.956 | 0.477 | 0.70 | 0.370 | 0.317 | 14.32
8467.254 | 0.567 | 0.44 | 0.439 | 0.376 | 14.35
8502.483 | 0.639 | 1.22 | 0.494 | 0.450 | 8.91
8665.019 | 1.188 | 1.83 | 0.913 | 0.839 | 8.11
8750.472 | 1.518 | 1.62 | 1.163 | 1.068 | 8.17
8862.782 | 1.925 | 2.99 | 1.468 | 1.389 | 5.38
Table 6: SMC N77A Hydrogen Spectrum De-Reddened and Residuals Wavelength | Obs. Flux | Error | De–Red | De–Red |  $|$ Diff. $|$
---|---|---|---|---|---
(Å) | (H$\beta$= 100) | 1$\sigma$ % | Obs. Flux | Model | %
3734.375 | 2.445 | 2.05 | 2.558 | 2.401 | 6.14
3750.158 | 2.950 | 1.32 | 3.084 | 3.059 | 0.81
3770.637 | 3.739 | 1.02 | 3.906 | 3.980 | 1.89
3835.391 | 6.730 | 0.66 | 7.012 | 7.336 | 4.62
3970.079 | 14.78 | 0.62 | 15.32 | 15.97 | 4.24
4101.742 | 25.04 | 0.34 | 25.82 | 25.99 | 0.66
4340.471 | 43.30 | 0.37 | 44.23 | 47.00 | 6.26
4861.333 | 100.00 | 0.59 | 100.00 | 100.00 | 0.00
6562.819 | 369.47 | 0.26 | 350.20 | 284.45 | 18.77
8437.956 | 0.334 | 2.91 | 0.307 | 0.307 | 0.15
8467.254 | 0.423 | 2.13 | 0.389 | 0.365 | 6.17
8502.483 | 0.418 | 3.72 | 0.385 | 0.437 | 13.51
8665.019 | 0.826 | 1.32 | 0.758 | 0.817 | 7.78
8750.472 | 1.135 | 1.98 | 1.041 | 1.041 | 0.03
8862.782 | 1.547 | 0.94 | 1.418 | 1.354 | 4.51
Table 7: SMC N88A Hydrogen Spectrum De-Reddened and Residuals Wavelength | Obs. Flux | Error | De–Red | De–Red |  $|$ Diff. $|$
---|---|---|---|---|---
(Å) | (H$\beta$= 100) | 1$\sigma$ % | Obs. Flux | Model | %
3734.375 | 2.010 | 1.04 | 2.591 | 2.392 | 7.68
3750.158 | 2.531 | 0.79 | 3.253 | 3.050 | 6.24
3770.637 | 2.994 | 0.78 | 3.832 | 3.970 | 3.60
3835.391 | 5.696 | 0.59 | 7.196 | 7.324 | 1.78
3970.079 | 13.20 | 0.32 | 16.23 | 15.95 | 1.73
4101.742 | 22.89 | 0.16 | 27.38 | 25.97 | 5.15
4340.471 | 42.41 | 0.27 | 48.23 | 47.13 | 2.28
4861.333 | 100.00 | 0.36 | 100.00 | 100.00 | 0.00
6562.819 | 471.06 | 0.21 | 327.41 | 283.07 | 13.54
8437.956 | 0.727 | 0.70 | 0.316 | 0.300 | 5.06
8467.254 | 0.881 | 0.26 | 0.381 | 0.355 | 6.82
8502.483 | 1.040 | 0.70 | 0.446 | 0.426 | 4.48
8665.019 | 1.832 | 0.25 | 0.759 | 0.795 | 4.74
8750.472 | 2.398 | 0.23 | 0.977 | 1.012 | 3.58
8862.782 | 3.147 | 0.13 | 1.255 | 1.316 | 4.86
Table 8: SMC N456 Hydrogen Spectrum De-Reddened and Residuals Wavelength | Obs. Flux | Error | De–Red | De–Red |  $|$ Diff. $|$
---|---|---|---|---|---
(Å) | (H$\beta$= 100) | 1$\sigma$ % | Obs. Flux | Model | %
3734.375 | 2.334 | 0.54 | 2.556 | 2.401 | 6.06
3750.158 | 2.989 | 0.84 | 3.269 | 3.059 | 6.42
3770.637 | 3.812 | 0.62 | 4.162 | 3.981 | 4.35
3835.391 | 6.726 | 0.47 | 7.307 | 7.339 | 0.44
3970.079 | 14.54 | 0.24 | 15.64 | 15.98 | 2.17
4101.742 | 24.20 | 0.13 | 25.76 | 26.00 | 0.93
4340.471 | 43.36 | 0.44 | 45.32 | 47.05 | 3.82
4861.333 | 100.00 | 0.43 | 100.00 | 100.00 | 0.00
6562.819 | 341.59 | 0.43 | 304.33 | 283.79 | 6.75
8437.956 | 0.439 | 21.00 | 0.356 | 0.305 | 14.33
8467.254 | 0.542 | 1.11 | 0.438 | 0.362 | 17.35
8502.483 | 0.577 | 13.50 | 0.466 | 0.434 | 6.87
8665.019 | 0.976 | 4.73 | 0.783 | 0.811 | 3.58
8750.472 | 1.403 | 0.51 | 1.123 | 1.032 | 8.10
8862.782 | 1.343 | 2.36 | 1.070 | 1.344 | 25.61
Table 9: Global Spherical Isobaric Derived Reddening Parameters. | Rv | $C_{n}$
---|---|---
Hii Region |  | (dex)
LMC N191 | 3.41 | 0.200
SMC N77A | 2.74 | 0.061
SMC N88A | 6.49 | 0.771
SMC N456 | 4.02 | 0.174
Table 10: Global Spherical Isobaric Nebula Derived Properties | $\log(Q_{H})$ | $T_{\rm src.}$ | $\log(P/k)$ | Fe Depl. | $<T_{\rm e}>$ | $<n_{\rm e}>$
---|---|---|---|---|---|---
Hii Region | (cm s-1) | (K) |  | (dex) | (K) | (cm-3)
LMC N191 | 7.625 | 37 500 | 6.80 | 0.5 | 9 274 | 339
SMC N77A | 8.125 | 38 500 | 6.00 | 1.0 | 11 860 | 43
SMC N88A | 9.250 | 46 500 | 7.40 | 1.5 | 14 360 | 1 016
SMC N456 | 8.125 | 38 500 | 6.45 | 1.5 | 12 540 | 114

####  6.3.1 Best-fit model

Table 9 presents the $R_{v}$ and $C_{n}$ for each nebula. The $R_{v}$ of SMC
N77A and LMC N191 are consistent with the $R_{v}=$2.93 of the Small Magellanic
Cloud and $R_{v}=$3.16 of the Large Magellanic Cloud (Pei, 1992) within an
error of $\sim$6%. The $R_{v}$ of SMC N88A and SMC N456 are larger than the
$R_{v}$ of the SMC.

Table 10 presents the derived parameters of the best-fit model. These four
nebulae have the central stellar temperature of 37 500 K$<T_{\rm eff}<$46 500
K. The ionization parameter of the nebulae in SMC is 8.125$<log(Q_{\rm
H})<$9.250, which are larger than the typical ionization parameter
7$<log(Q_{\rm H}<$8 of local Hii regions (Dopita et al., 2000). The nebula in
LMC has an ionization parameter of 7.625 falling in the range of ionization
parameters of typical local Hii regions. The ISM pressure is
6.0$<log(P/k)<$7.4 of these four nebulae.

The electron temperature of nebula derived from the best-fit model is 9 000
K$<T_{e}<$14 500 K and consistent with the volume-weighted average
temperature, $T_{e}$[O iii], measured from the IFU data. The temperature of
LMC N191 is 9 270 K, matching with the measured temperature of 8 864$\pm$813
K. The temperatures of SMC N77A and SMC N456 are close to the typical
assumption of $T_{e}=$10 000 K of an ideal ”Case B” nebula model. The SMC N77A
has the temperature of $T_{e}=$11 860 K consistent with the measured
temperature $T_{e}=$11 405$\pm$630 K. The SMC N456 has the derived temperature
of $T_{e}=$12 540 K in agreement with the measured temperature $T_{e}=$11
3215$\pm$423 K. The SMC 88A has the hottest temperature in the sample, with a
derived temperature $T_{e}=$14 170 K consistent with the measured temperature
of $T_{e}=$12 932$\pm$1 763 K

The density of the nebula derived from the best-fit model is 10$\rm
cm^{-3}$$<n_{e}<$$10^{3}$$\rm cm^{-3}$. The SMC N77A has the lowest density of
$n_{e}=$43$\rm cm^{-3}$. The LMC N191 and SMC N456 have the density of
$n_{e}=$339$\rm cm^{-3}$ and 114$\rm cm^{-3}$ respectively. The SMC N88A is
the densest nebula with a density of $n_{e}=$1016$\rm cm^{-3}$.

###  6.4 Emission-line Ratios

Table 11: Key line-ratios of the nebulae and the best-fit models. | [O iii]/[O ii] | [O iii]/[O iii] | [O ii]/[O ii] | [S ii]/[S ii]
---|---|---|---|---
| (5007+4959)/(3726+3729) | 4363/(5007+4959) | 3726/3729 | 6716/6731
LMC-N191 | 0.332 | 0.0035 | 0.90 | 1.15
$\sigma$% | 1.19 | 0.99 | 0.92 | 0.22
Model | 0.291 | 0.0041 | 0.93 | 1.08
$\Delta\%$ | 12.3 | 17.1 | 3.3 | 6.1
SMC N77A | 1.378 | 0.00754 | 0.724 | 1.344
$\sigma$% | 0.97 | 0.81 | 0.75 | 0.42
Model | 1.179 | 0.00800 | 0.719 | 1.385
$\Delta\%$ | 14.45 | 6.13 | 0.71 | 3.07
SMC N88A | 15.144 | 0.01415 | 1.236 | 0.875
$\sigma$% | 1.86 | 0.34 | 1.84 | 0.5
Model | 13.499 | 0.0140 | 1.232 | 0.849
$\Delta\%$ | 10.86 | 0.71 | 0.32 | 2.97
SMC N456 | 1.899 | 0.00795 | 0.759 | 1.298
$\sigma$% | 0.93 | 0.45 | 0.86 | 0.19
Model | 1.604 | 0.00987 | 0.767 | 1.301
$\Delta\%$ | 15.55 | 24.04 | 1.13 | 0.2

Table 11 presents the comparison of the key integrated line-ratios between the
observations and the best-fit model. We compare the ionization-parameter-
sensitive line-ratio [O iii]$\lambda\lambda$ 4959,5007/[O
ii]$\lambda\lambda$3726,27 (the well-known $O32$ ratio), the temperature-
sensitive line-ratio [O iii] $\lambda$4363/[O iii]$\lambda\lambda$4959,5007,
and the density-sensitive line-ratios [O ii]$\lambda$3726/[O ii]$\lambda$3729
and [S ii]$\lambda$6716 /[S ii]$\lambda$6731\. We use $\Delta\%$ to measure
the derivation of the best-fit line-ratios from the observation. $\Delta\%$ is
described as below,

$\Delta\%=\frac{\lvert
F_{\lambda}^{Obs}-F_{\lambda}^{Mod}\rvert}{F_{\lambda}^{Obs}}\times 100,$ (6)

where $F_{\lambda}^{Obs}$ and $F_{\lambda}^{Mod}$ are observed and modeled
fluxes.

The density-sensitive [O ii]-ratio and [S ii]-ratio are well-fitted with the
$\Delta\%$ smaller than 10%. The goodness of the temperature-sensitive line-
ratio fitting varies from the $\Delta\%$ of 0.71% in SMC N88A to the
$\Delta\%$ of 24.04% in SMC N456. The $\Delta\%$ of the [O iii]/[O ii]ratio
ranges from 10-15%.

Figure 10 presents the comparison of the radial profile of the temperature-
sensitive line-ratio, [O iii]$\lambda$4363/[O iii] $\lambda$5007, and the
density-sensitive line-ratios, [O ii]$\lambda 3727$ /[O ii]$\lambda$3729 and
[S ii]$\lambda 6717$/[S ii]$\lambda$6731, between the best-fit model and the
observation. The distance to the center of nebula is normalized by the maximal
radius of the nebula, $R_{max}$. We project the 3D spherical nebula onto a 2D
plane in order to compare with the observation. The line-ratio at each
distance is the integration along the line-of-sight direction.

The constant pressure model predicts that the [O ii]$\lambda 3727$/[O
ii]$\lambda$3729 ratio is an increasing function of the distance to the center
of nebula. In contrast, the observation shows that the [O ii]-ratio decreases
along the distance to the center of nebula for LMC N191 and SMC N456. In SMC
N77A and SMC N88A, the median value of the [O ii]-ratio displays a decreasing
trend along the distance to the nebular center although the measurement error
of 25-50% on the [O ii]-ratio is too large to reveal an evident gradient of
the radial profile.

The observations show that the [S ii]$\lambda 6717$ /[S ii]$\lambda$6731 ratio
increases with the distance to the center of nebula for the four Hii regions.
In contrast, the modeled [S ii]-ratio displays a flat distribution within the
nebula. The discrepancy of the [S ii]-ratio between the model and the
observation reveals that the density structure of the nebula is more complex
than the current models.

The temperature-sensitive [O iii]$\lambda$4363/[O iii] $\lambda$5007 ratio
profiles are flat in the four Hii regions, which are consistent with the
measurements of the flat temperature gradients in nebula. In contrast, the
constant pressure models predict a decreasing [O iii]-ratio profile for the
four Hii regions.

###  6.5 The pressure of nebula

Figure 11 presents the radial profile of the pressure of the four Hii regions.
For each nebula, the pressure, $log(P/k)$, is calculated by using four
different indicators.

  * •

$\log(P/k)^{OF06}_{\rm[SII]}$ is calculated by using $T_{e}$ measured by the
Osterbrock & Ferland (2006) method and $n_{e\rm[SII]}$ derived from [S
ii]-ratio by using the Osterbrock & Ferland (2006) method,

  * •

$\log(P/k)^{K19}_{\rm[SII]}$ is calculated by using $T_{e}$ measured by the
Osterbrock & Ferland (2006) method and $n_{e\rm[SII]}$ derived from [S
ii]-ratio by using the Kewley et al. (2019a) method,

  * •

$\log(P/k)^{OF06}_{\rm[OII]}$ is calculated by using $T_{e}$ measured by the
Osterbrock & Ferland (2006) method and $n_{e\rm[OII]}$ derived from [O
ii]-ratio by the Osterbrock & Ferland (2006) method,

  * •

$\log(P/k)^{K19}_{\rm[OII]}$ is calculated by using $T_{e}$ measured by the
Osterbrock & Ferland (2006) method and $n_{e\rm[OII]}$ derived from [O
ii]-ratio by the Kewley et al. (2019a) method,

The radial profiles of $log(P/k)$ shows an isobaric condition within the
nebula. In SMC N77A, SMC N88A and SMC N456, the pressure has a flat
distribution across the nebula. In LMC N191, the pressure has a flat gradient
within 0.7$R_{max}$ and decreases from 0.7$R_{max}$ to 1$R_{max}$. The average
pressure of the nebula ranges from 5.4$\pm$1.1$<log(P/k)<$7.7$\pm$0.36, which
is consistent with the range of pressures predicted by the best-fit models.

Figure 10: The comparison of radial distribution of the density-sensitive and
the temperature-sensitive line ratios between our observational data and our
constant pressure models. Each row presents the radial distribution of each
individual Hii region. The line ratios are the projection along the line-of-
sight directions. The radius is normalized to the maximum radius of models for
the constant pressure models and is normalized to the maximum radius of the
nebula for the observed data. Left: The comparison of the radial distribution
of the density-sensitive line ratio of [S ii]$\lambda$6716/[S
ii]$\lambda$6731\. Middle: The comparison of the radial distribution of the
density-sensitive line ratio of [O ii]$\lambda$3727/[O ii]$\lambda$3729\.
Left: The comparison of the radial distribution of the temperature-sensitive
line ratio of [O iii]$\lambda$5007/[O iii]$\lambda$4363\. In each panel, the
black points are the observed data and the lines present the radial
distribution of [S ii]$\lambda$6716/[S ii]$\lambda$6731 given by constant
pressure models with different log($P/k$) and ionization parameters.  Figure
11: The ISM pressure profiles of the Hii regions. The pressure is derived
through $P/k=2n_{e}(1+He/H)T_{e}$. The radius is normalized to the maximum
radius of models for the constant pressure models and is normalized to the
maximum radius of the nebula for the observed data. The red dashed line
indicates the average of the derived pressure. The purple dashed line
indicates the modeled pressure value. The yellow region indicates the
1$\sigma$ of the measured pressure.

##  7 Discussion

###  7.1 Temperature Structures of Hii regions

Figure 12: Comparison of the electron temperature of high-ionization species
and the temperature of low-ionization species. The temperature of high-
ionization species is measured from the [O iii]$\lambda$5007/[O
iii]$\lambda$4363 ratio. The temperature of low-ionization species is measured
from the [N ii]$\lambda\lambda$6548,84/[N ii]$\lambda$5755 ratio. The solid
line is the relationship given by López-Sánchez et al. (2012) and the dot-
dashed line is the relationship induced by Garnett (1992).

The electron temperature is a crucial diagnostic for determining the ISM
metallicity. Calibrations of metallicity are derived from simple
photoionization models with constant temperature or density. However, real Hii
regions are composed of multi-ionization zones with fluctuations of
temperatures and densities (Garnett, 1992). In our data, Hii regions have two
distinct zones: a high-ionization zone and a low-ionization zone. The high-
ionization zone with high-ionization species, like $\rm Ne^{2+}$ and $\rm
O^{2+}$, is located closer to the central star than the low-ionization zone
with low-ionization species, like $\rm O^{+}$ and $\rm N^{+}$, within nebula.
This is because the stellar radiation field is absorbed in regions closer to
the central star so that there are few ionizing photons to ionize atoms to
high ionization stages in the outer regions of the nebula..

Analyses of integrated nebular spectra suggest that electron temperatures vary
from a low-ionization zone to a high-ionization zone (Peimbert, 1967; Peimbert
et al., 2004, 2017). The lines [O iii]$\lambda 4363$, [O iii]$\lambda\lambda
4959,5007$ are most often used for determining the temperature in Hii regions.
However, the temperature derived from the [O iii] lines, $T_{e}\rm(OIII)$,
only represents the temperature in the $\rm O^{2+}$ zones. Hägele et al.
(2008) suggest that the missing $T_{e}$ in low-ionization zones can lead to an
underestimate of the metallicity by 0.2 dex. A two-ionization zone model is
assumed to convert the $T_{e}$ in high-ionization zones to the $T_{e}$ in low
ionization zones. As shown in Figure 12, our Hii region sample suggests that
the observed relationship of average temperatures in high-ionization zones and
low-ionization zones matches with the analytical function given by previous
research (Garnett, 1992; López-Sánchez et al., 2012).

In contrast to the assumption of constant temperature in photoionization
models, real Hii regions have radial variations of electron temperatures.
Detailed study of the $\theta^{1}$ Ori C area in the Orion nebula demonstrates
that the electron temperature rises with distance from the ionizing star
(Rubin et al., 2011). Our data show a sign of the positive gradient of the
temperature profiles within LMC N191, SMC N77A and SMC N88A. However, the
uncertainties of the best-fit gradients are too large to exclude the
possibility of a flat or negative gradient of the temperature profiles within
these Hii regions. SMC N456 shows a sign of a negative gradient of the
temperature profile but the uncertainty of the best-fit gradient is too large
to exclude a positive or flat gradient. The electron temperature gradient
causes an underestimate of the derived metallicity calibration (Stasińska,
2005) and can cause the difference between metallicities calibrated from
recombination lines and auroral lines (Peimbert et al., 2017).

The inhomogeneity of the electron temperature within Hii regions is attributed
to several main mechanisms (Peimbert et al., 2017). One mechanism is the
extreme density inhomogeneity leading to the low-temperature heavily shadowed
area within Hii regions. The second mechanism is the temperature inhomogeneity
is the mechanical energy input by stellar winds and shocks (Arthur et al.,
2016; O’Dell et al., 2017). Another potential cause of temperature
inhomogeneities is the spatial distribution of multiple ionizing sources
within Hii regions.

###  7.2 Density Structures of Hii regions

Most Hii region models assume constant density across the nebula. However, the
constant density assumption is not realistic given that real Hii regions have
significant density variations. The density variations are found to be 80-700
$\rm cm^{-3}$ in the Orion nebula (Rubin et al., 2011) and to be 40-4000 $\rm
cm^{-3}$ within Hii regions in the inner region of the Milky Way (Simpson et
al., 2004).

The density radial gradients are complex in nearby Hii regions. The
measurement of electron densities across the $\theta^{1}$ Ori C area shows
that the density variation is not a monotonic function of distance to the
ionizing star (Rubin et al., 2011). In large nearby Hii regions, the density
gradients are observed to be flat (García-Benito et al., 2010; Ramos-Larios et
al., 2010) while the density gradients appear to be negative in compact Hii
regions (Binette et al., 2002; McLeod et al., 2016). Particularly, ultra-
compact Hii regions have the steepest density gradients among all categories
of nebulae (Kurtz, 2002; Johnson & Kobulnicky, 2003; Phillips, 2007).

Our observed Hii regions also show a large range of density gradients. The
density gradient is as steep as around $-700~{}\rm cm^{-3}~{}pc^{-1}$ in SMC
N88A. However, the density gradients can also be as flat as $-3~{}\rm cm^{-3}$
$\rm pc^{-1}$ in SMC N456.

Diverse density gradients of Hii regions are attributed to complex nebular
geometry. In the Orion nebula, complex density gradients are likely the result
of a turbulent Hii region (Arthur et al., 2016; Ha et al., 2021) with bar-like
structures (van der Werf et al., 2013; Rubin et al., 2011), where the
contamination of scattered light from dense clumps change the electron density
distribution (Kewley et al., 2019b).

There are two major processes proposed that create complex density structures
within Hii regions. One is the “collect and collapse” scenario (Elmegreen et
al., 1995), where dense clumps are formed when Hii regions expand into
turbulent ISM. The other is the “radiation driven implosion” model (Bertoldi,
1989), where overdensities within Hii regions are amplified by the heating of
hot stars. Both processes create an overdensity of the ISM (Walch et al.,
2015; Schneider et al., 2016). The stellar wind is another potential cause for
the complex geometries of Hii regions (Park et al., 2010).

###  7.3 Pressure structures of Hii regions

The ISM pressure is a key parameter to describe the ISM properties, which
includes both nebular temperature and density structures. The ISM pressure is
determined by the mechanical energy produced by stellar feedback, the strength
and the shape of the radiation field. Models with a constant ISM pressure are
more realistic than models with a constant density or temperature, when the
sound-crossing timescale is shorter than the heating and cooling timescale.
Previous research shows that the condition of constant pressure occurs within
the majority of Hii regions (Begelman & McKee, 1990; Gutiérrez & Beckman,
2010).

The value of ISM pressure crosses a broad interval in Hii region models
(Kewley et al., 2019a, b), ranging from a low pressure of log($P/k$)=4 to a
high pressure environment with log($P/k$)=8-9 which is observed in high-
redshift Hii regions (Lehnert et al., 2009). Our sampled Hii regions in the
LMC and the SMC present ISM pressures around log($P/k$)=6-8, which is
consistent with the range given by models.

###  7.4 Implications for future photoionization models

Self-consistent photoionization models, like cloudy and mappings, require as
input an assumption of the shape of the nebula. Constant temperature or
density assumptions are always applied to simple photoionization models but
this is not physically realistic. Our observations reveal that the
substructures, as filaments or dense knots, exist within real Hii regions,
which are caused by stellar winds, ISM turbulence and the nebular expansion.
Our observations indicate that the temperature and density structures do not
obey the isothermal or constant density assumptions. Gradients of temperature
and density structures exist in the LMC and SMC nebulae. Allowing temperature
and density to vary results in a more physically realistic model.

To model realistic Hii regions, some approximations of nebular geometries are
applied. Some studies have used plane-parallel approximations to mimic the
behavior of non-spherical structures in nebula (Levesque et al., 2010), and
some use concentric shells to deal with the shell-like structures within Hii
regions (Pellegrini et al., 2020). In these models, the contribution of the
diffuse ionizing photons to the ionization field is assumed as a fixed
fraction, which may bear no relation to a real nebula with complex geometries.

Fully self-consistent three-dimensional photoionization codes are needed.
Monte-Carlo radiative transfer techniques offer the promise of a substantial
improvement over current simple structures to handle the complex geometry of
nebulae. Some Monte-Carlo photoionization codes are already available
(Ercolano & Bastian, 2008; Vandenbroucke & Wood, 2018) and have shown promise
in modeling the ionized gas around young stars (Law et al., 2011) and nebulae
with filamentary structures (Ercolano et al., 2012). For a precise modeling of
nebular emission-lines and internal structures, a Monte-Carlo code with
comprehensive considerations of atomic data and ISM microphysics, like cooling
and heating processes, is under development. The new code incorporates the
Monte-Carlo radiative transfer technique into the existing mappings
photoionization code and will be presented in a forthcoming paper (Jin et al.
in prep.)

This research was supported by the Australian Research Council Centre of
Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D), through
project number CE170100013. L.J.K. gratefully acknowledges the support of an
ARC Laureate Fellowship (FL150100113).

## References

  * Arellano-Córdova & Rodríguez (2020) Arellano-Córdova, K. Z. & Rodríguez, M. 2020, MNRAS, 497, 672. doi:10.1093/mnras/staa1759
  * Arthur et al. (2016) Arthur, S. J., Medina, S.-N. X., & Henney, W. J. 2016, MNRAS, 463, 2864
  * Begelman & McKee (1990) Begelman, M. C. & McKee, C. F. 1990, ApJ, 358, 375. doi:10.1086/168994
  * Bertoldi (1989) Bertoldi, F. 1989, ApJ, 346, 735. doi:10.1086/168055
  * Binette (1985) Binette, L. 1985, A&A, 143, 334
  * Binette et al. (2002) Binette, L., González-Gómez, D. I., & Mayya, Y. D. 2002, Rev. Mexicana Astron. Astrofis., 38, 279
  * Childress et al. (2014) Childress, M. J., Vogt, F. P. A., Nielsen, J., et al. 2014, Ap&SS, 349, 617
  * Cioni (2009) Cioni, M.-R. L. 2009, A&A, 506, 1137
  * Chu & Kennicutt (1986) Chu, Y.-H., & Kennicutt, R. C. 1986, ApJ, 311, 85
  * Copetti et al. (2000) Copetti, M. V. F., Mallmann, J. A. H., Schmidt, A. A., et al. 2000, A&A, 357, 621
  * Del Zanna et al. (2015) Del Zanna, G., Dere, K. P., Young, P. R., et al. 2015, A&A, 582, A56
  * Dopita et al. (2000) Dopita, M. A., Kewley, L. J., Heisler, C. A., et al. 2000, ApJ, 542, 224. doi:10.1086/309538
  * Dopita et al. (2007) Dopita, M., Hart, J., McGregor, P., et al. 2007, Ap&SS, 310, 255
  * Dopita et al. (2010) Dopita, M., Rhee, J., Farage, C., et al. 2010, Ap&SS, 327, 245
  * Dopita et al. (2013) Dopita, M. A., Sutherland, R. S., Nicholls, D. C., et al. 2013, ApJS, 208, 10
  * Dopita et al. (2016) Dopita, M. A., Seitenzahl, I. R., Sutherland, R. S., et al. 2016, ApJ, 826, 150. doi:10.3847/0004-637X/826/2/150
  * Dopita et al. (2018) Dopita, M. A., Vogt, F. P. A., Sutherland, R. S., et al. 2018, ApJS, 237, 10. doi:10.3847/1538-4365/aac837
  * Dopita & Seitenzahl (2019) Dopita, M. A., & Seitenzahl, I. 2019, Supernova Remnants: An Odyssey in Space After Stellar Death II, 4
  * Elmegreen et al. (1995) Elmegreen, B. G., Kimura, T., & Tosa, M. 1995, ApJ, 451, 675. doi:10.1086/176253
  * Ercolano & Bastian (2008) Ercolano, B., & Bastian, N. 2008, Mass Loss from Stars and the Evolution of Stellar Clusters, 397
  * Ercolano et al. (2012) Ercolano, B., Dale, J. E., Gritschneder, M., et al. 2012, MNRAS, 420, 141
  * Field (1965) Field, G. B. 1965, ApJ, 142, 531. doi:10.1086/148317
  * Fitzpatrick (1999) Fitzpatrick, E. L. 1999, PASP, 111, 63
  * Filippenko (1985) Filippenko, A. V. 1985, ApJ, 289, 475
  * García-Benito et al. (2010) García-Benito, R., Díaz, A., Hägele, G. F., et al. 2010, MNRAS, 408, 2234
  * Garnett (1992) Garnett, D. R. 1992, AJ, 103, 1330. doi:10.1086/116146
  * Garnett et al. (1995) Garnett, D. R., Skillman, E. D., Dufour, R. J., et al. 1995, ApJ, 443, 64. doi:10.1086/175503
  * Garnett & Dinerstein (2001) Garnett, D. R., & Dinerstein, H. L. 2001, ApJ, 558, 145
  * Gordon et al. (2003) Gordon, K. D., Clayton, G. C., Misselt, K. A., et al. 2003, ApJ, 594, 279. doi:10.1086/376774
  * Groves et al. (2004) Groves, B. A., Dopita, M. A., & Sutherland, R. S. 2004, ApJS, 153, 75
  * Gutiérrez & Beckman (2010) Gutiérrez, L. & Beckman, J. E. 2010, ApJ, 710, L44. doi:10.1088/2041-8205/710/1/L44
  * Ha et al. (2021) Ha, T., Li, Y., Xu, S., et al. 2021, ApJ, 907, L40. doi:10.3847/2041-8213/abd8c9
  * Hägele et al. (2006) Hägele, G. F., Pérez-Montero, E., Díaz, Á. I., et al. 2006, MNRAS, 372, 293
  * Hägele et al. (2008) Hägele, G. F., Díaz, Á. I., Terlevich, E., et al. 2008, MNRAS, 383, 209
  * Henize (1956) Henize, K. G. 1956, ApJS, 2, 315. doi:10.1086/190025
  * Herrera-Camus et al. (2016) Herrera-Camus, R., Bolatto, A., Smith, J. D., et al. 2016, ApJ, 826, 175
  * Hubeny & Lanz (1995) Hubeny, I., & Lanz, T. 1995, ApJ, 439, 875
  * Hubeny & Lanz (2017) Hubeny, I., & Lanz, T. 2017, arXiv e-prints, arXiv:1706.01859
  * Jin et al. (2022) Jin, Y., Kewley, L. J., & Sutherland, R. S. 2022, ApJ, 934, L8. doi:10.3847/2041-8213/ac80f3
  * Johnson & Kobulnicky (2003) Johnson, K. E. & Kobulnicky, H. A. 2003, ApJ, 597, 923. doi:10.1086/378585
  * Johnstone et al. (2015) Johnstone, C. P., Güdel, M., Lüftinger, T., et al. 2015, A&A, 577, A27. doi:10.1051/0004-6361/20142530010.48550/arXiv.1503.06669
  * Kennicutt (1984) Kennicutt, R. C. 1984, ApJ, 287, 116
  * Kennicutt & Hodge (1986) Kennicutt, R. C., & Hodge, P. W. 1986, ApJ, 306, 130
  * Kewley et al. (2019a) Kewley, L. J., Nicholls, D. C., Sutherland, R., et al. 2019, ApJ, 880, 16
  * Kewley et al. (2019b) Kewley, L. J., Nicholls, D. C., & Sutherland, R. S. 2019, ARA&A, 57, 511
  * Kollmeier et al. (2017) Kollmeier, J. A., Zasowski, G., Rix, H.-W., et al. 2017, arXiv e-prints, arXiv:1711.03234
  * Krabbe & Copetti (2002) Krabbe, A. C., & Copetti, M. V. F. 2002, A&A, 387, 295
  * Kumari et al. (2017) Kumari, N., James, B. L., & Irwin, M. J. 2017, MNRAS, 470, 4618
  * Kurtz (2002) Kurtz, S. 2002, Hot Star Workshop III: The Earliest Phases of Massive Star Birth, 267, 81
  * Law et al. (2011) Law, K.-H., Gordon, K. D., & Misselt, K. A. 2011, ApJ, 738, 124
  * Lehnert et al. (2009) Lehnert, M. D., Nesvadba, N. P. H., Le Tiran, L., et al. 2009, ApJ, 699, 1660
  * Levesque et al. (2010) Levesque, E. M., Kewley, L. J., & Larson, K. L. 2010, AJ, 139, 712
  * Liu et al. (2006) Liu, X.-W., Barlow, M. J., Zhang, Y., et al. 2006, MNRAS, 368, 1959
  * López-Sánchez et al. (2012) López-Sánchez, Á. R., Dopita, M. A., Kewley, L. J., et al. 2012, MNRAS, 426, 2630. doi:10.1111/j.1365-2966.2012.21145.x
  * Malmann et al. (2002) Malmann, J. A. H., Copetti, M. V. F., Dos Santos, I. C. F., et al. 2002, Revista Mexicana De Astronomia Y Astrofisica Conference Series, 57
  * Mendoza (1983) Mendoza, C. 1983, Planetary Nebulae, 103, 143
  * McLeod et al. (2016) McLeod, A. F., Gritschneder, M., Dale, J. E., et al. 2016, MNRAS, 462, 3537. doi:10.1093/mnras/stw1864
  * Nicholls et al. (2020) Nicholls, D. C., Kewley, L. J., & Sutherland, R. S. 2020, PASP, 132, 033001. doi:10.1088/1538-3873/ab6818
  * O’Dell et al. (2017) O’Dell, C. R., Ferland, G. J., & Peimbert, M. 2017, MNRAS, 464, 4835
  * Osterbrock (1989) Osterbrock, D. E. 1989, Astrophysics of Gaseous Nebulae and Active Galactic Nuclei
  * Osterbrock & Ferland (2006) Osterbrock, D. E., & Ferland, G. J. 2006, Astrophysics of gaseous nebulae and active galactic nuclei
  * Park et al. (2010) Park, J.-W., Min, K.-W., Seon, K.-I., et al. 2010, ApJ, 719, 1964. doi:10.1088/0004-637X/719/2/1964
  * Pei (1992) Pei, Y. C. 1992, ApJ, 395, 130. doi:10.1086/171637
  * Peimbert (1967) Peimbert, M. 1967, ApJ, 150, 825
  * Peimbert & Torres-Peimbert (1976) Peimbert, M. & Torres-Peimbert, S. 1976, ApJ, 203, 581. doi:10.1086/154114
  * Peimbert et al. (2004) Peimbert, M., Peimbert, A., Ruiz, M. T., et al. 2004, ApJS, 150, 431
  * Peimbert et al. (2017) Peimbert, M., Peimbert, A., & Delgado-Inglada, G. 2017, PASP, 129, 082001. doi:10.1088/1538-3873/aa72c3
  * Pellegrini et al. (2020) Pellegrini, E. W., Rahner, D., Reissl, S., et al. 2020, MNRAS, doi:10.1093/mnras/staa1473
  * Peña-Guerrero et al. (2012a) Peña-Guerrero, M. A., Peimbert, A., Peimbert, M., et al. 2012, ApJ, 746, 115
  * Peña-Guerrero et al. (2012b) Peña-Guerrero, M. A., Peimbert, A., & Peimbert, M. 2012, ApJ, 756, L14
  * Phillips (2007) Phillips, J. P. 2007, MNRAS, 380, 369. doi:10.1111/j.1365-2966.2007.12078.x
  * Ramos-Larios et al. (2010) Ramos-Larios, G., Phillips, J. P., & Pérez-Grana, J. A. 2010, MNRAS, 405, 245. doi:10.1111/j.1365-2966.2010.16483.x
  * Rubin et al. (1994) Rubin, R. H., Simpson, J. P., Lord, S. D., et al. 1994, ApJ, 420, 772
  * Rubin et al. (2011) Rubin, R. H., Simpson, J. P., O’Dell, C. R., et al. 2011, MNRAS, 410, 1320. doi:10.1111/j.1365-2966.2010.17522.x
  * Rubin et al. (2016) Rubin, R. H., Simpson, J. P., Colgan, S. W. J., et al. 2016, MNRAS, 459, 1875
  * Russell & Dopita (1992) Russell, S. C., & Dopita, M. A. 1992, ApJ, 384, 508
  * Schneider et al. (2016) Schneider, N., Bontemps, S., Motte, F., et al. 2016, A&A, 591, A40. doi:10.1051/0004-6361/201628328
  * Selier & Heydari-Malayeri (2012) Selier, R., & Heydari-Malayeri, M. 2012, A&A, 545, A29
  * Simpson et al. (2004) Simpson, J. P., Rubin, R. H., Colgan, S. W. J., et al. 2004, ApJ, 611, 338. doi:10.1086/422028
  * Stasińska (2005) Stasińska, G. 2005, A&A, 434, 507
  * Stasińska et al. (2013) Stasińska, G., Morisset, C., Simón-Díaz, S., et al. 2013, A&A, 551, A82
  * Sutherland & Dopita (1993) Sutherland, R. S., & Dopita, M. A. 1993, ApJS, 88, 253
  * Testor & Lortet (1987) Testor, G. & Lortet, M.-C. 1987, A&A, 178, 25
  * Testor et al. (2003) Testor, G., Lemaire, J. L., & Field, D. 2003, A&A, 407, 905
  * Testor et al. (2010) Testor, G., Lemaire, J. L., Heydari-Malayeri, M., et al. 2010, A&A, 510, A95. doi:10.1051/0004-6361/200911786
  * Toribio San Cipriano et al. (2017) Toribio San Cipriano, L., Domínguez-Guzmán, G., Esteban, C., et al. 2017, MNRAS, 467, 3759. doi:10.1093/mnras/stx328
  * van der Werf et al. (2013) van der Werf, P. P., Goss, W. M., & O’Dell, C. R. 2013, ApJ, 762, 101. doi:10.1088/0004-637X/762/2/101
  * Vandenbroucke & Wood (2018) Vandenbroucke, B., & Wood, K. 2018, Astronomy and Computing, 23, 40
  * Vermeij & van der Hulst (2002) Vermeij, R., & van der Hulst, J. M. 2002, A&A, 391, 1081
  * Walch et al. (2015) Walch, S., Girichidis, P., Naab, T., et al. 2015, MNRAS, 454, 238. doi:10.1093/mnras/stv1975