text:
\section{Introduction}
\label{sec:intro}

\begin{table}
	\begin{center}
		\begin{tabular}{|c|l|}
			\hline
		 \textbf{Type of Operations}   & \textbf{Parameters} \\ \hline \hline
		Median Filtering       & kernel: 7x7, 5x5, 3x3     \\  \hline
	    Gaussian Blurring      & std: 1.1, 0.75, 0.5       \\ \hline
	    Additive Noise  (AWGN) & std: 2.0, 0.5, 0.25       \\ \hline
		Resizing               & factor: 1.5, 1.125, 1.01  \\ \hline
		JPEG Compression       & quality: 1.5, 1.125, 1.01  \\ \hline
		\end{tabular}
	\end{center}
	\caption{List of considered image manipulations}
	\label{tab:type}
\end{table}

 problema della forgery detection.
Stato dell'arte..

\begin{figure}
	\centering
	\includegraphics[width=0.9\linewidth]{figure/spamLD}
	\caption{Classic scheme of Rich model features}
	\label{schemaFridrichLD}
\end{figure}

\begin{figure}
	\centering
	\includegraphics[width=0.9\linewidth]{figure/LD}
	\caption{Classic scheme of Local Description}
	\label{schemaLD}
\end{figure}

\begin{figure}
	\centering
	\includegraphics[width=0.9\linewidth]{figure/cnnLD}
	\caption{Classic scheme of Local Description}
	\label{schemaCnnLD}
\end{figure}

\begin{table}
	\begin{center}
		\begin{tabular}{|c|c||c||c|c|}
			\cline{3-5}
          \multicolumn{2}{c|}{} & Bayer2016 & SPAM & CNN   \\
          \multicolumn{2}{c|}{} & 60 epoch & +SVM & 15 epoch \\ \hline \hline
			Median      &  7x7  &        &  99.61 &         \\ \cline{2-5}
			Filtering   &  5x5  &  96.66 &  99.67 &  99.47  \\ \cline{2-5}
			            &  3x3  &  94.56 &  99.83 &  99.35  \\ \hline \hline
			Gaussian    &  1.1  &  99.65 &  99.93 &  99.79  \\ \cline{2-5}
			Blurring    &  0.75 &        &  99.90 &         \\ \cline{2-5}
			            &  0.5  &  83.10 &  87.10 &  95.70  \\ \hline \hline
			Additive    &  2.0  &  97.08 &  99.94 &  99.95  \\ \cline{2-5}
			Noise       &  0.5  &        &  99.37 &         \\ \cline{2-5}
			            &  0.25 &  51.83 &  85.06 &  88.81  \\ \hline \hline
			Resizing    &  1.5  &  99.22 &  99.99 & 100.00  \\ \cline{2-5}
			            & 1.125 &  91.06 &  98.94 &  99.56  \\ \cline{2-5}
			            &  1.01 &  80.51 &  96.01 &  97.81  \\ \hline \hline
			JPEG        &   70  &        &  99.99 &         \\ \cline{2-5}
			Compression &   80  &  77.01 &  99.73 &  99.37  \\ \cline{2-5}
			            &   90  &  63.77 &  90.86 &  92.08  \\ \hline
		\end{tabular}
	\end{center}
	\caption{Detection Accuracy for binary classification}
	\label{tab:perf}
\end{table}

\begin{table}
	\begin{center}
		\begin{tabular}{|c||c|c|c|c|c|}
			\cline{2-6}
	\multicolumn{1}{c|}{}  & Med. & Blur & Noise & Res. & JPEG \\ \hline \hline
	Median Filtering       & 99.30 &  0.40 &  0.09 &  0.12 &  0.09  \\  \hline
	Gaussian Blurring      &  1.71 & 90.96 &  0.17 &  6.83 &  0.34  \\ \hline
	Additive Noise         &  0.15 &  0.44 & 98.04 &  0.68 &  0.69  \\ \hline
	Resizing               &  0.25 &  7.77 &  0.01 & 90.64 &  1.33  \\ \hline
	JPEG Compression       &  0.06 &  0.50 &  0.51 &  0.85 & 98.09  \\ \hline
		\end{tabular}
	\end{center}
	\caption{Confusion matrix SPAM+SVM}
	\label{tab:conf_spam}
\end{table}

\section{Conclusions}

\bibliographystyle{IEEEbib}

\section{Introduction}
\label{sec:intro}

Images and videos represent by now the dominant source of traffic on the Internet and the bulk of data stored on social media.
Nowadays, however, such data may be easily manipulated by malicious attackers
to convey some twisted and potentially dangerous messages.
Interested areas include politics, journalism, judiciary, even the scientific world.
For this reason, in the last few years there has been intense and ever growing activity on multimedia forensics, aiming at developing methods to detect, localize, and classify possible image
manipulations.

Typical attacks consist in adding or deleting objects,
using material taken from the same image (copy-move) or from other sources (splicing).
As a consequence, many researchers have focused on detecting near-duplicates in the image or across a repository of images.
More fundamentally, one may be interested in establishing
whether the image under analysis is pristine or else has been subjected to some post-processing,
since any form of manipulation may raise suspects and suggest deeper inquiry.
In fact, most of the times,
copy-moves and splicing are accompanied by various forms of elaboration aimed at removing the most
obvious traces of editing.
These include, for example, resizing, rotation, linear and non-linear filtering,
contrast enhancement, histogram equalization and, eventually, re-compression.

A number of papers have been proposed to detect one or the other of such elaborations
\cite{Popescu2005, Kirchner2008, Kirchner2010, Huang2010, Stamm2010}.
These methods, however, are sensitive to just some specific manipulations.
A more appealing line of research is to detect {\em all} possible manipulations,
an approach that has been followed in several papers \cite{Cao2012, Verdoliva2014, Fan2015, Li2016}.
Notably, in the 2013 IEEE Image Forensics Challenge,
the most effective techniques for both image forgery detection \cite{Cozzolino2014a} and localization \cite{Cozzolino2014b}
used this approach, relying on powerful residual-based local descriptors.
These features,
such as SPAM (subtractive pixel adjacency matrix) \cite{Pevny2010} or SRM (spatial rich models) \cite{Fridrich2012},
inspired to previous work in steganalysis, are extracted from the so-called residual image.
In fact, the noise residual, extracted through some high-pass filtering of the image,
contains a wealth of information on the in-camera and out-camera processes involved in the image formation.
Such subtle traces, hardly visible without enhancement,
may reveal anomalies due to object insertion \cite{Cozzolino2015b, Cozzolino2016}
or can detect different types of image editing operations \cite{Verdoliva2014, Li2016, Boroumand2017}.

Very recently,
inspired by impressive results in the closely related fields of computer vision and pattern recognition \cite{Krizhevsky2012},
the multimedia forensics community began focusing on the use of deep learning \cite{Bengio2013, LeCun2015},
especially convolutional neural networks (CNN) \cite{Zeiler2014}.
Taking advantage of the lesson learnt from SPAM/SRM features,
constrained CNN architectures have been proposed both for steganalysis \cite{Qian2015} and manipulation detection \cite{Bayar2016},
where the first convolutional layer is forced to perform a high-pass filtering.

In this paper we show that there is no real contraposition between residual-based features and CNNs.
Indeed, these local features can be computed through a CNN with architecture and parameters selected so as to guarantee a perfect equivalence.
Once established this result,
we go beyond emulation, removing constraints on parameters, and fine-tuning the net to further improve performance.
Since the resulting network has a lightweight structure,
fine-tuning can be carried out by means of a limited training set, limiting computation time and memory usage.
A significant performance gain with respect to the conventional feature is observed, especially in the most challenging situations.

In the following we describe in more detail the residual-based local features (Section 2),
recast them as a constrained CNN, to be further trained after removing constraints (Section 3),
show experimental results (Section 4),
and draw conclusions (Section 5).

\section{Residual-based local descriptors}

Establishing which type of processing an image has undergone
calls for the ability to detect the subtle traces left by these operations,
typically in the form of recurrent micropatterns.
This problem has close ties with steganalysis, where weak messages hidden in the data are sought,
so it is no surprise that the same tools, residual-based local descriptors,
prove successful in both cases.
To associate a residual-based feature to an image, or an image block, the following processing chain
has been successfully used in steganalysis \cite{Pevny2010, Fridrich2012}:
\begin{enumerate}
	\item   extraction of noise residuals
	\item   scalar quantization
	\item   computation of co-occurrences
	\item   computation of histogram
\end{enumerate}

In the following, we describe in some more depth all these steps,
taking a specific model out of the 39 proposed in \cite{Fridrich2012} as running example.

\vspace{3mm} \noindent
{\bf Extraction of noise residual.}
The goal is to extract image details, in the high-frequency part of the image,
which enables the analysis of expressive micropatterns.
As the name suggests, this step can be implemented by resorting to a high-pass filter.
In \cite{Fridrich2012} a number of different high-pass filters have been considered,
both linear and nonlinear, with various supports.
Here, as an example,
we focus on a single 4-tap mono-dimensional linear filter, with coefficients
${\bf w} = [1, -3, 3, -1]$.
The filter extracts image details along one direction, but is applied also on the image transpose (assuming vertical/horizontal invariance) to augment the available data.
Choosing a single filter rather than considering all models proposed in \cite{Fridrich2012}
is motivated not only by the reduced complexity
but also by the very good performance observed in the context of image forgery detection \cite{Cozzolino2014a, Verdoliva2014, Li2016}.

\vspace{3mm} \noindent
{\bf Scalar quantization.}
Residuals are conceptually real-valued quantities or, in any case,
high-resolution integers, so they must be quantized to reduce cardinality and allow easy processing.
In \cite{Fridrich2012} a uniform quantization is used with an odd number of levels (to ensure that 0 is among the possible outputs).
Therefore, the only parameters to set are the number of quantization levels, $L$, and the quantization step $\Delta$.
In our example we set $L=3$ and $\Delta=4.5$.

\vspace{3mm} \noindent
{\bf Co-occurrences.}
The computation of co-occurrences is the core step of the procedure.
In fact, this is a low-complexity means for taking into account high-order dependencies among residuals
and hence gather information on recurrent micropatterns.
Following \cite{Fridrich2012} we compute co-occurrences on $N=4$ pixels in a row,
both along and across the filter direction.
With these values, two co-occurrence matrices with $3^4=81$ entries are obtained.
Of course, the image or block under analysis must be large enough to obtain meaningful estimates.
All the co-occurrence $N$-dimensional bins are eventually coded as integers.

\vspace{3mm} \noindent
{\bf Feature formation.}
Counting co-occurrences one obtains the final feature vector describing the image.
Neglecting symmetries, our final feature has length equal to 162.
The final classification phase is performed by a linear SVM classifier.

\section{Recasting local features as CNN}

We will now show that local residual-based features can be extracted by means of a convolutional neural network.
Establishing this equivalence leaves us with a CNN architecture and a set of parameters that are already known to provide an excellent performance for the problem of interest.
Then, given this good starting point,
we can move a step forward and fine-tune the network through a sensible training phase with labeled data.
Note that in this way we will carry out a joint optimization of both the feature extraction process and classification.
In the following we will first move from local features to a Bag-of-Words (BoW) paradigm,
and then proceed to the implementation by means of Convolutional Neural Networks.

\begin{figure}[t!]
    \centering
    \includegraphics[width=1.0\linewidth]{figure/schema00to01.pdf}
    \caption{Basic processing scheme (0) for extracting the single model SRM feature,
    and equivalent scheme (1) with inverted order of scalar quantization (SQ) and $n$-pixel shifting ($z^{-n}$).}
    \label{schemaFridrichLD}
\end{figure}

\subsection{From local features to Bag-of-Words}

In Fig.1 (top) we show the basic processing scheme used to extract the single model SRM feature.
Let $\X$ be the input image\footnote{We use capital boldface for images, lowercase boldface for vectors, and simple lowercase for scalars.
The value of image $\X$ at spatial site $s$, will be denoted by $x_s$.}, $\R$ the residual image, and $\hR$ the quantized residual image.
To compute the output feature, the input image is high-pass filtered,
then the residual image is quantized, and $N$ versions of it are generated, shifted one pixel apart from one another.
For each pixel $s$, the values $\hr_{1,s},\ldots,\hr_{N,s}$ are regarded as base-$L$ digits and encoded as a single scalar $i^*_s$,
finally, the histogram ${\bf h}$ of this latter image is computed.
The scheme at the bottom of Fig.1 is identical to the former except for the inverted order of scalar quantization (SQ) and shifting.
This inversion, however, allows us to focus on the two groups of blocks highlighted at the top of Fig.2.

The filter-shifter group can be replaced by a bank of $N$ filters, all identical to one another except for the position of the non-zero weights.
So, with reference to our running example,
the $n$-th filter will have non-zero weights, [1, -3, 3, -1], only on the $n$-th row, and zero weights everywhere else.
Turning to the second group,
the combination of $N$ scalar quantizers can be regarded as a constrained form of vector quantization (VQ).
More specifically, it is a product VQ, since the VQ codebook is obtained as the cartesian product of the $N$ SQ codebooks.
On one hand, product quantization is much simpler and faster than VQ.
On the other hand, its strong constraints are potentially detrimental for performance.
Its $K=L^N$ codewords are forced to lie on a truncated $N$-dimensional square lattice \cite{Gray1998} and cannot adapt to the data distribution.
Many of them will be wasted in empty regions of the feature space,
causing a sure, and often severe, loss of performance with respect to unconstrained VQ.

However, the most interesting observation about the new structure at the bottom of Fig.2
is that it implements the Bag-of-Words (or also Bag-of-Features) paradigm.
The filter bank extracts a feature vector for each image pixel, based on its neighborhood.
These feature are then associated, through VQ, with some template features.
Finally, the frequency of occurrence of the latter, computed in the last block, provides a synthetic descriptor of the input image.
The fact that filters and vector quantizer are largely sub-optimal impacts only on performance, not on interpretation.
Needless to say, they could be both improved through supervised training.

\begin{figure}[t!]
	\centering
	\includegraphics[width=1.0\linewidth]{figure/schema01to02}
	\caption{The cascade of filter and shifters of scheme (1) can be replaced by a bank of filters,
    while the bank of independent SQ's + coder can be replaced by product VQ.
    The resulting scheme (2) fits the Bag-of-Words paradigm.}
	\label{schemaFridrichLD}
\end{figure}

\begin{figure}[t]
	\centering
	\includegraphics[width=0.9\linewidth]{figure/schemaVQ2}
	\caption{A vector quantizer (left) can be implemented through a bank of filter followed by argmax (right).}
	\label{schemaLD}
\end{figure}

\subsection{From Bag-of-Words to CNN}

We now show that the processing steps of Fig.2 can be all implemented through a CNN.
First of all, the bank of linear filters used to extract noise residuals can be replaced by a pure convolutional layer,
with neurons computing the residuals as
\begin{equation}
    r_{n,s} = f(w_{n,s} * x_s + b_n)   \hspace{5mm} n=1,\cdots,N
\end{equation}
with $s$ used for image spatial location and $n$ to identify neurons.
The neuron weights coincide with filter coefficients, biases $b_n$ are all set to zero, and the non-linearity $f(\cdot)$ is set to identity.

As for the vector quantizer,
assuming the usual minimum distance hard-decision rule, it can be implemented by means of a convolutional layer followed by an hard-max layer.
Let $\r_s$ be a vector formed by collecting a group of residuals at site $s$, and $\ci_k$ the $k$-th codeword of the quantizer.
Their squared Euclidean distance, $d^2_{k,s}$, can be expanded as
\begin{eqnarray}
    ||\r_s-\ci_k||^2  & = & ||\r_s||^2 + ||\ci_k||^2  - 2<\r_s,\ci_k> \nonumber\\
                      & = & ||\r_s||^2 - 2 \epsilon_k - 2<\r_s,\ci_k>
\end{eqnarray}
with $||\cdot||^2$ and $<\cdot,\cdot>$ indicating norm and inner product, respectively.
Hence, neglecting the irrelevant $||\r_s||^2$ term:
\begin{eqnarray}
    i^*_s & = & {\rm argmin}_{k=1,\cdots,K} \,\, d_{k,s} \nonumber\\
          & = & {\rm argmax}_{k=1,\cdots,K} \,\,(<\r_s,\ci_k> + \epsilon_k) \nonumber\\
          & = & {\rm argmax}_{k=1,\cdots,K} \,\, m_{k,s}
\end{eqnarray}
with $m_{k,s}$ interpreted as a matching score between the feature vector at site $s$ and the $k$-th codeword.

\begin{figure}[t!]
	\centering
	\includegraphics[width=1\linewidth]{figure/schema02to03}
	\caption{The whole scheme (2) can be converted in the CNN (3).
    The filter bank is replaced by a convolutional layer, VQ is replaced by convolutional-hardmax layers, the histogram can be computed through an average pooling layer.}
	\label{schemaLD}
\end{figure}

This equivalence is depicted in Fig.3.
The matching scores $m_{k,s}$ are computed through a convolutional layer,
equipped with $K=L^N$ filters (remember that $L$ is the number of quantization levels),
one for each codeword, having weights $\ci_k$, bias $\bm{\epsilon}_k$ and,
again, an identity as activation function.
The best matching codeword is then selected through a hardmax processing.

The evolution of the whole network is shown in Fig.4.
As already said, the first filter bank is replaced by a convolutional layer.
Then, the VQ is replaced by another convolutional layer followed by a hardmax layer.
The former outputs $K$ feature maps, ${\bf M}_k$, with the matching scores.
The latter outputs $K$ more binary maps, ${\bf P}_k$,
where $p_{k,s}=1$ when the corresponding matching score $m_{k,s}$ is maximum over $k$, and 0 otherwise.
Finally, the histogram computation is replaced by an average pooling layer operating on the whole feature maps, that is, $h_k = \sum_s p_{k,s}$.
The resulting scheme is shown at the bottom of Fig.4.

\begin{figure}
	\centering
	\includegraphics[width=1\linewidth]{figure/schema03to04}
	\caption{The constrained CNN of scheme (3) extracts the features which feed an external classifier.
    In scheme (4) this is replaced by an internal fully connected layer,
    and all constraints are removed.
    By fine-tuning on training data, all layers can be optimized jointly.
    The dashed lines are to remind that this net provides only half the feature, a twin net (not shown for clarity) provides the other half.}
	\label{schemaCnnLD}
\end{figure}

This net, actually, computes only half the desired feature, the part based on across-filter co-occurrences.
The twin half is computed similarly, and the complete feature is eventually fed to the SVM classifier, as shown in Fig.5 (top).
In this network, weights and biases of the convolutional layers
are all hard-wired to reproduce exactly the behavior of the residual-based local feature described in Section 2,
thereby ensuring the good performance observed in the literature.
We now proceed to remove all constraints and allow the net to learn on a suitable training set.
First of all, the classifier itself can be implemented as part of the CNN architecture by including a fully connected layer at the end,
obtaining the architecture of Fig.5 (bottom).
Now, to exploit the full potential of deep learning, all parameters must be optimized by appropriate training,
thereby overcoming all the impairing constraints mentioned before.
Note that the learning phase allows us not only to optimize all layers, which could be done also in the BoW framework,
but to optimize them {\em jointly}, taking full advantage of the CNN structural freedom.
Moreover, the lightweight architecture of the network
is instrumental to achieve good results even in the presence of a limited training set.

ayer with a sft-max layer that approximates it, so as to avoid non-differentiable operators.

However, before proceeding with the training, it is necessary to replace the hard-max layer,
with a soft-max layer that approximates it, so as to avoid non-differentiable operators.
Given the input vector $\{m_{k,s}, k=1,\ldots,K\}$, the soft-max computes the quantities
\begin{equation}
    p_{k,s} = e^{\alpha m_{k,s}}/\sum_l e^{\alpha m_{l,s}}
\end{equation}
With the aim to preserve a close correspondence with the original descriptor, we should choose a very large $\alpha$ parameter, so as to obtain a steep nonlinearity.
However, as said before, this is not really necessary, since our goal is only to improve performance.
Hence, we select a relatively small value for $\alpha$ in order not to slow down learning.
Likewise, to implement minimum distance VQ exactly, the biases in the second convolutional layer should depend on the filter weights,
but there is no practical reason to enforce this constraint, and we allow also the biases to adapt freely.

Now, the final CNN can be trained as usual with stochastic gradient descent \cite{Krizhevsky2012} to adapt to the desired task.
It should be clear that a number of architectural modifications could be also tested starting from this basic structure,
but this goes beyond the scope of the present paper, and will be the object of future research.
We must underline that the equivalence between CNN and BoW has been noticed before in the literature,
for example in \cite{Lan2015, Arandjelovic2016, Richard2016}.

\section{Experimental analysis}

\renewcommand{\tabcolsep}{6pt}
\begin{table}
    \caption{Image manipulations under test.}
	\begin{center}
		\begin{tabular}{|l|l|} \hline
			\ru			Manipulation           & Parameters                                 \\ \hline \hline
			\ru			Median Filtering       & kernel: 7$\times$7, 5$\times$5, 3$\times$3 \\ \hline
			\ru			Gaussian Blurring      & st. dev: 1.1,\; 0.75,\; 0.5                \\ \hline
			\ru			Additive Noise  (AWGN) & st. dev.: 2.0,\; 0.5,\; 0.25               \\ \hline
			\ru			Resizing               & scale: 1.5,\; 1.125,\; 1.01                \\ \hline
			\ru			JPEG Compression       & quality factor: 70,\;  80,\; 90            \\ \hline
		\end{tabular}
	\end{center}
	\label{tab:type}
\end{table}

\renewcommand{\tabcolsep}{6pt} \begin{table*}
	\label{tab:perf}
	\caption{Detection Accuracy for binary classification tasks.}
	\begin{center}
		\begin{tabular}{|c|c||c|c|c||c|c|c|}  \hline
			\multicolumn{2}{|c||}{\ru Manipulation} & \multicolumn{3}{c||}{Small Training Set} & \multicolumn{3}{c|}{Large Training Set} \\ \hline
			\multicolumn{2}{|c||}{\ru}              &   Bayar2016 &     SRM+SVM &    prop. CNN &   Bayar2016 &     SRM+SVM &   prop. CNN \\
			\multicolumn{2}{|c||}{   }              &   60 epochs &             &    15 epochs &   60 epochs &             &   15 epochs \\ \hline \hline
			                  & \ru    7x7          &       98.23 &       99.61 &       99.07  &       99.69 &       99.68 &       99.55 \\ \cline{2-8}
			Median Filtering  & \ru    5x5          &       96.66 &       99.67 &       99.47  &       99.78 &       99.68 &       99.60 \\ \cline{2-8}
			                  & \ru    3x3          &       94.56 &       99.83 &       99.35  &       99.80 &       99.87 &       99.75 \\ \hline \hline
			                  & \ru    1.1          &       99.65 &       99.93 &       99.79  &       99.98 &       99.97 &       99.95 \\ \cline{2-8}
			Gaussian Blurring & \ru    0.75         &       98.52 &       99.90 &       99.82  &       99.94 &       99.77 &       99.93 \\ \cline{2-8}
			                  & \ru    0.5          &       83.10 &       87.10 &       95.70  &       94.57 &       87.55 &       96.56 \\ \hline \hline
			                  & \ru    2.0          &       97.08 &       99.94 &       99.95  &       99.56 &       99.94 &       99.94 \\ \cline{2-8}
			Additive Noise    & \ru    0.5          &       82.93 &       99.37 &       99.36  &       93.83 &       99.34 &       99.66 \\ \cline{2-8}
			                  & \ru    0.25         &       51.83 &       85.06 &       88.81  &       80.28 &       84.01 &       90.79 \\ \hline \hline
			                  & \ru    1.5          &       99.22 &       99.99 &      100.00  &       99.72 &       99.87 &      100.00 \\ \cline{2-8}
			Resizing          & \ru ~~1.125~~       &       91.06 &       98.94 &       99.56  &       97.02 &       96.00 &       99.78 \\ \cline{2-8}
			                  & \ru    1.01         &       80.51 &       96.01 &       97.81  &       98.44 &       95.11 &       97.20 \\ \hline \hline
			                  & \ru     70          &       96.04 &       99.99 &       99.99  &       99.43 &       99.99 &       99.94 \\ \cline{2-8}
			JPEG Compression  & \ru     80          &       77.01 &       99.73 &       99.37  &       98.12 &       99.94 &       99.86 \\ \cline{2-8}
			                  & \ru     90          &       63.77 &       90.86 &       92.08  &       79.69 &       90.90 &       94.59  \\ \hline
		\end{tabular}
	\end{center}
\end{table*}

To test the performance of the proposed CNN architecture we carry out a number of experiments with typical manipulations.
Our synthetic dataset includes images taken from 9 devices,
4 smartphones (Apple iPhone 4S, Apple iPhone 5s, Huawei P7 mini, Nokia Lumia 925) and 5 cameras (Canon EOS 450D, Canon IXUS 95 IS, Sony DSC-S780, Samsung Digimax 301, Nikon Coolpix S5100).
Each device contributes 200 images, and from each image non-overlapping patches of dimension $128\times 128$ are sampled.
We select at random 6 devices to form the training set, while the remaining 3 are used as a testing set.
Therefore, the patches used for testing come from devices that are never seen in the training phase.
For each pristine patch, the corresponding manipulated patch is also included in the set.
Overall, our training set comprises a total of just 10800+10800 patches, quite a small number for deep learning applications.
We consider 5 types of image manipulation: median filtering, gaussian blurring, AWGN noise addition, resizing, and JPEG compression,
with three different settings for each case (see Tab.1) corresponding to increasingly challenging tasks.
For example, JPEG compression with quality factor Q=70 is always easily detected,
while a quality factor Q=90 makes things much more difficult.

In the proposed CNN
the first convolutional layer includes 4 filters of size $5\times5\times1$ operating on the monochrome input (we use only the green band normalized in $[0,1]$.
In the second layer there are 81 filters of size $1\times1\times4$.
Filters are initialized as described in Section 3 and the $\alpha$ parameter of soft-max is set to $2^{16}$.
The code is implemented in Tensorfow and runs on a Nvidia Tesla P100 with 16GB RAM.
We set the learning rate to $10^{-6}$, with decay $5 \cdot 10^{-4}$, batch size 36 and Adam \cite{Kingma2015} optimization method, using the cross-entropy loss function.
Together with the proposed CNN we consider also the basic solution,
with the handcrafted feature followed by linear SVM,
and the CNN proposed in Bayar2016 \cite{Bayar2016}
based on the use of a preliminary high-pass convolutional layer.

Results in terms of probability of correct decision for each binary classification problem
are reported in the left part of Tab.2 (small training set).
With ``easy'' manipulations, {\it e.g.} JPEG@70, all methods provide near-perfect results
and there is no point in replacing the SRM+SVM solution with something else.
In the presence of more challenging attacks, however, the performance varies significantly across methods.
After just 15 epochs of fine tuning, the proposed CNN improves over SRM+SVM of about 2 percent points for JPEG compression, resizing, and noising, and more than 8 points for blurring,
while median filtering is almost always detected in any case.
In the same cases, the CNN architecture proposed in \cite{Bayar2016} provides worse results,
sometimes close to 50\%, even after 60 epochs of training.
Our conjecture is that a deep CNN is simply not able to adapt correctly with a small training set.
In this condition a good hand-crafted feature can work much better.
The proposed CNN builds upon this result and takes advantage of the available limited training data
to fine-tune its parameters.

To carry out a fair comparison we also considered a case in which a much larger training set is available,
comprising 460800 patches, that is more than 20 times larger than before.
Results are reported in the right part of Tab.2.
As expected the performance does not change much for the SRM+SVM solution,
since the SVM needs limited training anyway.
For the proposed CNN, some improvements are observed for the more challenging tasks.
As an example, for JPEG@90 the accuracy grows from 92.08 to 94.59.
Much larger improvements are observed for the network proposed in  \cite{Bayar2016},
which closes almost always the performance gap and sometimes outperforms slightly the proposed CNN.
Nonetheless in a few challenging cases,
like the already mentioned JPEG@90 or the addition of low-power white noise,
there is still a difference of more than 10 percent point with our proposal.
It is also interesting to compare the two adjacent columns,
proposed CNN at 15 epochs and the CNN architecture proposed in \cite{Bayar2016} at 60 epochs,
which speak clearly in favor of the first solution, in terms of both complexity and performance.

\begin{figure*}[t!]
	\centering
	\begin{tabular}{cccc}
		\includegraphics[width=0.22\linewidth]{figure/img/test3_ori.jpg} \;
		\includegraphics[width=0.22\linewidth]{figure/img/test3_forg.png} \;
		\includegraphics[width=0.22\linewidth]{figure/img/test3_050_gaussian05_SPAM.png} \;
		\includegraphics[width=0.22\linewidth]{figure/img/test3_050_gaussian05.png}  \\ \vspace{-3mm} \\
		\includegraphics[width=0.22\linewidth]{figure/img/test2_ori.jpg} \;
		\includegraphics[width=0.22\linewidth]{figure/img/test2_forg.png} \;
		\includegraphics[width=0.22\linewidth]{figure/img/test2_113_resize1125_SPAM.png} \;
		\includegraphics[width=0.22\linewidth]{figure/img/test2_113_resize1125.png}
	\end{tabular}
	\caption{Top: splicing (blurred with st. dev. 0.5). Bottom: copy-move (resized with scale 1.125).
		From left to right: original image, forged image, SRM+SVM heat map, proposed CNN heat map.}
	\label{fig:localization1}
\end{figure*}

The ability to reliably classify small
patches may be very valuable in the presence of spatially localized attacks.
This is the case of image copy-move or splicing, where only a small part of the image is tampered with.
In these cases, a descriptor computed on small patches can more reliably detect manipulations,
and even localize the forgery by working in sliding-window modality.
Fig.6 shows two examples of forgery localization with a slightly blurred splicing and a resized copy-move, respectively.
In both cases
using the proposed CNN in sliding-window modality, a sharp heat map is obtained,
allowing for a precise localization of the forgery.
The SRM+SVM solution also provides good results, but the heat map is more fuzzy, with a higher risk of false alarms.
It is worth underlining again that the images used for these tests come from cameras that did not contribute to the training set.

\section{Conclusions}

Residual-based descriptors have proven extremely effective for a number of image forensic applications.
Improving upon the current state of the art, however, is slow and costly,
since the design of better hand-crafted features is not trivial.
We showed that a class of residual-based features can be regarded as compact constrained CNNs.
This represent a precious starting point to exploit the huge potential of deep learning,
as testified by the promising early results.
However, this is only a first step,
and there is much room for improvements, especially through new architectural solutions.
This will be the main focus of future work.

\balance

\vspace{3mm}
{\footnotesize
{\bf Acknowledgment.}
This material is based on research \vspace{-0.6mm}
sponsored by the Air Force Research Laboratory \vspace{-0.6mm}
and the Defense Advanced Research Projects Agency \vspace{-0.6mm}
under agreement number FA8750-16-1-0104.
The U.S. Government is authorized to reproduce and distribute reprints \vspace{-0.6mm}
for Governmental purposes notwithstanding any copyright notation thereon. \vspace{-0.6mm}
The views and conclusions contained herein are those of the authors \vspace{-0.6mm}
and should not be interpreted as necessarily representing the official policies or endorsements, \vspace{-0.6mm}
either expressed or implied, of the Air Force Research Laboratory \vspace{-0.6mm}
and the Defense Advanced Research Projects Agency or the U.S. Government.
}

\bibliographystyle{IEEEtran}

redpajama_set_name: RedPajamaArXiv


text:
\section{Introduction}
\label{sec:intro}

Argonium (ArH$^+$), the first noble gas molecule discovered in space, was detected by \citet{barlow2013} in emission in the J=1-0 and J=2-1 transitions in the Crab Nebula, along with an OH$^+$ transition at 971 GHz. The Crab Nebula is a young supernova remnant (SNR), consisting of a central pulsar wind nebula (PWN) and a surrounding network of filaments ionised by the PWN synchrotron emission \citep{hester2008}. Molecular hydrogen (H$_2$) \citep{graham1990,loh2010,loh2011,loh2012} emission had previously been detected in the Crab Nebula, with modelling by \citet{richardson2013} suggesting it originates from a trace molecular component in mostly atomic gas, shielded from ionising radiation at the centre of the filaments. \citet{barlow2013} noted that ArH$^+$ is principally formed by the reaction Ar$^+$ + H$_2$ $\to$ ArH$^+$ + H \citep{roach1970}, while the main destruction mechanism is ArH$^+$ + H$_2$ $\to$ Ar + H$_3^+$ \citep{schilke2014}.

Since its discovery in the Crab Nebula, ArH$^+$ has been found to be a fairly common molecule in the interstellar medium (ISM). \citet{schilke2014} reported ArH$^+$ J=1-0 absorption at a range of velocities along several galactic sightlines, while \citet{muller2015} detected ArH$^+$ absorption from a foreground galaxy in the spectrum of a lensed blazar. \citet{schilke2014} found the interstellar ArH$^+$ abundance to be uncorrelated with the abundances of the molecular gas tracers HF and H$_2$O$^+$, while chemical modelling of ArH$^+$ in the ISM by \citet{schilke2014} and \citet{neufeld2016} showed that it traces gas where the hydrogen is almost entirely atomic. ArH$^+$ was predicted to peak in abundance at the edges of interstellar clouds, where the ultraviolet (UV) radiation field is strongest.

The Crab Nebula is a substantially different environment to the interstellar clouds considered by \citet{schilke2014}. The PWN is a strong source of synchrotron radiation, including X-rays, which behave very differently to the UV radiation in the interstellar case. Additionally, the filaments must be permeated by the high energy charged particles which produce the synchrotron radiation, equivalent to but much more intense than the cosmic rays included in ISM astrochemical modelling. Both can provide an additional source of Ar$^+$, potentially allowing more ArH$^+$ to be formed, while the strong radiation field and generally harsh environment would tend to destroy molecules. This paper will attempt to model the chemistry of ArH$^+$ under the conditions present in the Crab Nebula with the goal of reproducing the observed emission.

\section{Description of \textsc{ucl\_pdr}}
\label{sec:code}

In this work we make use of the substantially upgraded 1D \textsc{ucl\_pdr} photodissociation region (PDR) code. It was first developed by \citet{bell2005,bell2006} and further augmented by \citet{bayet2011}. The code assumes a semi-infinite slab geometry, and the chemistry and thermal balance are calculated self-consistently at each depth point into the slab, producing chemical abundances, emission line strengths and gas temperatures as a function of depth. The gas temperature is calculated by taking into account a wide range of heating and cooling processes (see \citet{bell2006} and \citet{rollig2007}, where the code has also been benchmarked against other codes). The gas can be heated by the photoejection of electrons from dust grains and polycyclic aromatic hydrocarbons (PAHs), FUV pumping, photodissociation of H$_2$ molecules, cosmic rays, exothermic reactions, C I photoionization, H$_2$ formation, gas-grain collisions and turbulence. The gas is mainly cooled through emission from collisionally excited atoms and molecules and by interactions with the cooler dust grains. The original code included cooling from [O I], [C I] and [C II] fine structure lines and CO rotational lines (up to J=11--10). \citet{bayet2011} improved the \textsc{ucl\_pdr} code by adding the cooling, and hence the calculation of the line brightnesses, from $^{13}$CO, C$^{18}$O, CS and HCN transitions. The line brightnesses for all the coolant species are calculated using the Large Velocity Gradient (LVG) approximation (see e.g. \citet{vandertak2007}). All the collisional rates are taken from the Leiden Atomic and Molecular Database \citep{schoier2010}. The chemical network is user-specified, using the rate-coefficient formalism of the UMIST database \citep{mcelroy2013}.

The \textsc{ucl\_pdr} code has now been further improved by the following additions, made by T. Bell:
\begin{enumerate}
\item LVG calculation of line intensities (and cooling contributions) from arbitrary new species. This requires appropriate collisional data to be available.
\item Two-sided FUV illumination, either in the form of scaled Draine fields or full radiation SEDs which can be take from a radiative transfer model.
\item Photoreaction rates are now calculated from available cross-section data when full input radiation SEDs are provided.
\item A detailed treatment of the H$_2$ formation rate based on laboratory data, either in the form of a rate equation \citep{cazaux2004} or combined with results from advanced Monte Carlo rate calculations \citep{iqbal2012}.
\item The possibility to run gas-grain models including freeze-out, surface reactions, and various desorption mechanisms.
\item Implementation of the X-ray chemistry and physics, following \citet{meijerink2005} and \citet{stauber2005}. The shape and intensity of the incident X-ray spectrum can be specified, as appropriate for either YSOs or AGN.
\item Inclusion of a limited PAH chemistry, following \citet{wolfire2003,wolfire2008}.
\item Inclusion of an updated vibrationally excited H$_2$ chemistry and vibrational de-excitation heating rate calculation \citep{rollig2006}.
\item Fully parallel OpenMP implementation, with good scaling efficiency (the speed-up scales well with increasing number of cores).
\end{enumerate}
The user is now able to specify most model parameters within input text files, including dust extinction properties, input SEDs, cloud density structure, collisional data for LVG calculations, photorate cross sections, and so on. \textsc{ucl\_pdr} will shortly be made publically available from https://uclchem.github.io.

\section{Method}
\label{sec:method}

We have modelled a single Crab Nebula filament using a combination of two codes: \textsc{mocassin}, a Monte Carlo photoionisation code \citep{ercolano2003,ercolano2005,ercolano2008}, and the upgraded version of \textsc{ucl\_pdr} described above. \citet{owen2015} used \textsc{mocassin} to model the global properties of the Crab Nebula, fitting the predicted SED and line strengths to multi-wavelength observations in order to determine the dust and gas properties for a given distribution of matter. We used the dust and gas densities, dust grain size distribution and the elemental composition of the gas from their clumpy model VI, assuming amorphous carbon dust grains. We adopted a shell of matter with these properties located at 2.5 pc from a central point source for modelling with \textsc{mocassin}, using the Crab SED from \citet{hester2008} and the total SED luminosity of $1.3 \times 10^{38}$ erg s$^{-1}$ adopted by \citet{owen2015}. The shell thickness was adjusted so that the transition between ionised and neutral gas occurred at the outer edge, which required a thickness of $3.5 \times 10^{16}$ cm. The output SED was integrated between 912 \r{A} and 2400 \r{A} to give the transmitted UV field longward of the Lyman limit - $31$ Draines (see \citet{draine1978}) - and from 0.1 \r{A} to 100 \r{A} to give an X-ray flux of $0.35$ erg cm$^{-2}$ s$^{-1}$.

The filament was modelled in \textsc{ucl\_pdr} as a one dimensional slab, with the UV and X-ray fluxes that had been determined from \textsc{mocassin} incident on one side of the slab. The hydrogen nucleus number density, $n_{\rm{H}}$, was initially set to 1900 cm$^{-3}$, the value found for the clumps in model VI of \citet{owen2015}, while the initial cosmic ray ionisation rate was left as the standard interstellar value ($\zeta_0 = 1.3 \times 10^{-17} \rm{s}^{-1}$). The gas phase elemental abundances are listed in Table~\ref{tab:abun}. To account for the conditions in the Crab Nebula, \textsc{ucl\_pdr} was modified to include the effect of the Crab's higher dust-to-gas mass ratio ($\frac{M_d}{M_g} = 0.031$; \citet{owen2015}) on the H$_2$ formation and heating rates. The cosmic ray heating rate was taken from the code \textsc{cloudy} \citep{ferland1998}, as the default rate in \textsc{ucl\_pdr} (from \citet{goldsmith2001}) is for dark clouds, and may not be appropriate for the Crab Nebula, while the X-ray heating rate was calculated using the heating efficiency for pure helium, the most abundant element in the Crab, taken from \citet{dalgarno1999}. The X-ray absorption cross-sections were changed to those appropriate for the Crab Nebula abundances, and the input X-ray spectrum was approximated by a three component power law fit to the \textsc{mocassin} output spectrum, as shown in Figure~\ref{fig:xrayspec}.

\begin{table}
  \centering
  \caption{Gas phase elemental abundances, relative to hydrogen, used in PDR modelling, taken from model IV of \citet{owen2015}.}
  \begin{tabular}{cccc}
    \hline
    Element & Abundance & Element & Abundance \\
    \hline
    H & $1.00$ & N & $2.5 \times 10^{-4}$ \\
    He & $1.85$ & O & $6.2 \times 10^{-3}$ \\
    C & $1.02 \times 10^{-2}$ & Ar & $1.0 \times 10^{-5}$ \\
    \hline
  \end{tabular}
  \label{tab:abun}
\end{table}

\begin{figure}
  \centering
  \includegraphics[width=\columnwidth]{fig1.eps}
  \caption{The output X-ray spectrum from \textsc{mocassin} (solid line; see text) with the three component power law fit used by \textsc{ucl\_pdr} (dashed line).}
  \label{fig:xrayspec}
\end{figure}

ArH$^+$ was incorporated into an existing XDR chemical network using the reaction network from \citet{schilke2014}, with their photodissociation rate replaced by the one calculated by \citet{roueff2014} for the Crab Nebula, with the correction given by \citet{neufeld2016} (Ar + H$_2^+ \to$ ArH$^+$ + H, instead of Ar + H$_2^+ \to$ Ar$^+$ + H$_2$). Additional reactions to account for the effects of X-rays were included. The argon reactions incorporated and their rates are listed in Table~\ref{tab:xrayreac}. The rates of reactions 13 to 16 were calculated according to the methods detailed by \citet{meijerink2005}. The rates for reactions 17 and 18 were taken from \citet{shull1982} and from \citet{kingdon1996}, respectively. We also included reactions involving the H$^-$ anion from the UMIST database \citep{mcelroy2013}, to account for its role in the formation of H$_2$ in the Crab Nebula \citep{richardson2013}. Electron impact excitation rates for the first six rotational levels of ArH$^+$ were taken from \citet{hamilton2016}, and we used energy levels and Einstein A-coefficients from the Cologne Database for Molecular Spectroscopy (CDMS) \citep{muller2001,muller2005}, in order to calculate the emissivity of ArH$^+$ rotational transitions.

\begin{table*}
  \centering
  \caption{Argon reactions used for modelling the Crab nebula ArH$^+$ emission. CR refers to cosmic rays, CRPHOT to secondary photons produced by cosmic rays, XR to X-rays, XRSEC to secondary electrons produced by X-rays, XRPHOT to secondary photons from X-rays, h$\nu$ to a photon, $\zeta$ to cosmic ray or X-ray ionisation rates, $\sigma(E)$ is the X-ray ionisation cross section at energy $E$, $F(E)$ is the X-ray flux at energy $E$ and $\omega$ is the dust albedo.}
  \begin{tabular}{cccc}
    \hline
    & Reaction & Rate coefficient & Source \\
    \hline
    1 & Ar + CR $\to$ Ar$^{+}$ + e$^-$ & $10 \zeta_{\rm{H,cr}}$ & \citet{schilke2014} \\
    2 & Ar + CRPHOT $\to$ Ar$^+$ + e$^-$ & $0.8 \frac{\zeta_{\rm{H}_2\rm{,cr}}}{1 - \omega}$ & \citet{schilke2014} \\
    3 & Ar + H$_2^+$ $\to$ ArH$^+$ + H & $10^{-9}$ cm$^3$ s$^{-1}$ & \citet{schilke2014} \\
    4 & Ar + H$_3^+$ $\to$ ArH$^+$ + H$_2$ & $8 \times 10^{-10}$ exp$\left(\frac{6400 \rm{K}}{T}\right)$ cm$^3$ s$^{-1}$ & \citet{schilke2014} \\
    5 & Ar$^+$ + e$^-$ $\to$ Ar + h$\nu$ & $3.7 \times 10^{-12} \left(\frac{T}{300 \rm{K}}\right)^{-0.651}$ cm$^3$ s$^{-1}$ & \citet{schilke2014} \\
    6 & Ar$^+$ + H$_2$ $\to$ ArH$^+$ + H & $8.4 \times 10^{-10} \left(\frac{T}{300 \rm{K}}\right)^{0.16}$ cm$^3$ s$^{-1}$ & \citet{schilke2014} \\
    7 & ArH$^+$ + H$_2$ $\to$ Ar + H$_3^+$ & $8 \times 10^{-10}$ cm$^3$ s$^{-1}$ & \citet{schilke2014} \\
    8 & ArH$^+$ + CO $\to$ Ar + HCO$^+$ & $1.25 \times 10^{-9}$ cm$^3$ s$^{-1}$ & \citet{schilke2014} \\
    9 & ArH$^+$ + O $\to$ Ar + OH$^+$ & $8 \times 10^{-10}$ cm$^3$ s$^{-1}$ & \citet{schilke2014} \\
    10 & ArH$^+$ + C $\to$ Ar + CH$^+$ & $8 \times 10^{-10}$ cm$^3$ s$^{-1}$ & \citet{schilke2014} \\
    11 & ArH$^+$ + e$^-$ $\to$ Ar + H & $10^{-9}$ cm$^3$ s$^{-1}$ & \citet{schilke2014} \\
    12 & ArH$^+$ + h$\nu$ $\to$ Ar$^+$ + H & $4.20 \times 10^{-12}$ exp$\left(-3.27A_V\right)$ s$^{-1}$ & \citet{roueff2014} \\
    13 & Ar + XR $\to$ Ar$^{++}$ + e$^-$ + e$^-$ & $\int_{E_{\rm{min}}}^{E_{\rm{max}}} \sigma(E) \frac{F(E)}{E} \rm{d}E$ & \citet{verner1995} \\
    14 & Ar$^+$ + XR $\to$ Ar$^{++}$ + e$^-$ & $\int_{E_{\rm{min}}}^{E_{\rm{max}}} \sigma(E) \frac{F(E)}{E} \rm{d}E$ & \citet{verner1995} \\
    15 & Ar + XRSEC $\to$ Ar$^{+}$ + e$^-$ & $5.97 \zeta_{\rm{H,X-ray}}$ & \citet{lennon1988} \\
    16 & Ar + XRPHOT $\to$ Ar$^{+}$ + e$^-$ & $0.8 \frac{\zeta_{\rm{H}_2\rm{,X-ray}}}{1 - \omega}$ & \citet{schilke2014} \\
    17 & Ar$^{++}$ + e$^-$ $\to$ Ar$^+$ + h$\nu$ & $2.71 \times 10^{-11} \left(\frac{T}{300 \rm{K}}\right)^{-0.75}$ cm$^3$ s$^{-1}$ & \citet{shull1982} \\
    18 & Ar$^{++}$ + H $\to$ Ar$^+$ + H$^+$ & $10^{-15}$ cm$^3$ s$^{-1}$ & \citet{kingdon1996} \\
    \hline
  \end{tabular}
  \label{tab:xrayreac}
\end{table*}

\section{Results}
\label{sec:results}

As a check, \textsc{ucl\_pdr} was run using the same parameters as for the model interstellar cloud studied by \citet{schilke2014}. Figure~\ref{fig:schilke} shows the resulting abundances of Ar$^+$ and ArH$^+$ as a function of the visual extinction depth, $A_V$. The abundances of both species are found to be slightly smaller than those in \citet{schilke2014}, but are within a factor of a few throughout the cloud, and show the same behaviour with increasing $A_V$. As shown by \citet{schilke2014}, the main factor affecting the Ar$^+$ and ArH$^+$ abundances is the molecular hydrogen fraction, which reacts with and destroys both species.

\begin{figure}
  \centering
  \includegraphics[width=\columnwidth]{fig2.eps}
  \caption{Ar$^+$ (solid line) and ArH$^+$ (dashed line) abundances relative to hydrogen nuclei as a function of $A_V$, for an interstellar cloud model using the same parameters as \citet{schilke2014}.}
  \label{fig:schilke}
\end{figure}

Figure~\ref{fig:crababun} shows the abundances as a function of $A_V$ of H$_2$, Ar$^+$, ArH$^+$ and OH$^+$ (also observed by \citet{barlow2013}), for an initial Crab Nebula model clump. The model parameters used are listed in Table~\ref{tab:models}. The H$_2$ abundance is similar throughout the clump to the \citet{schilke2014} IS cloud model, but the abundance of ArH$^+$ is enhanced by over an order of magnitude, due to the additional production of Ar$^+$ by X-rays. The ArH$^+$ abundance is larger than or comparable to the peak interstellar cloud value in Figure~\ref{fig:schilke} for all but the highest $A_V$ values. The molecular ion OH$^+$, which has a similar formation mechanism (O$^+$, from cosmic ray or X-ray ionisation of neutral O, reacting with H$_2$), is also present at fairly high abundances, although its peak abundance is located at a higher $A_V$, and therefore deeper into the filament, than ArH$^+$.

\begin{table*}
  \centering
  \caption{Model parameters used for our initial Crab Nebula filament, and for our final models successfully reproducing the Herschel SPIRE FTS observations. $k$(e$^-$) is the rate coefficient of the ArH$^+$ + e$^-$ dissociative recombination reaction.}
  \begin{tabular}{cccccccc}
    \hline
    Model & $n_{\rm{H}}$ / cm$^{-3}$ & $\zeta$ / $\zeta_0$ & $F_{UV}$ / Draines & $F_X$ / erg cm$^{-2}$ s$^{-1}$ &  $A_V$ / mag & $A_V$/$N_{\rm{H}}$ / mag cm$^2$ & $k$(e$^-$) / cm$^3$ s$^{-1}$ \\
    \hline
    Initial & $1900$ & $1$ & $31$ & $0.35$ & $3.0$ & $6.289 \times 10^{-22}$ & $10^{-9}$ \\
    D3Z7 & $1900$ & $10^7$ & $31$ & $0.35$ & $0.05$ & $6.289 \times 10^{-22}$ & $10^{-11}$ \\
    D3Z8 & $1900$ & $10^8$ & $31$ & $0.35$ & $0.05$ & $6.289 \times 10^{-22}$ & $10^{-11}$ \\
    D4Z7 & $2 \times 10^4$ & $10^7$ & $31$ & $0.35$ & $0.01$ & $6.289 \times 10^{-22}$ & $10^{-11}$ \\
    D4Z8 & $2 \times 10^4$ & $10^8$ & $31$ & $0.35$ & $0.01$ & $6.289 \times 10^{-22}$ & $10^{-11}$ \\
    AVD3Z7 & $1900$ & $10^7$ & $31$ & $0.35$ & $0.30$ & $2.264 \times 10^{-21}$ & $10^{-11}$ \\
    AVD3Z8 & $1900$ & $10^8$ & $31$ & $0.35$ & $0.30$ & $2.264 \times 10^{-21}$ & $10^{-11}$ \\
    AVD4Z7 & $2 \times 10^4$ & $10^7$ & $31$ & $0.35$ & $0.30$ & $2.264 \times 10^{-21}$ & $10^{-11}$ \\
    AVD4Z8 & $2 \times 10^4$ & $10^8$ & $31$ & $0.35$ & $0.30$ & $2.264 \times 10^{-21}$ & $10^{-11}$ \\
    \hline
  \end{tabular}
  \label{tab:models}
\end{table*}

\begin{figure}
  \centering
  \includegraphics[width=\columnwidth]{fig3.eps}
  \caption{H$_2$, Ar$^+$, ArH$^+$ and OH$^+$ abundances as a function of $A_V$, for an initial Crab Nebula model. The H$_2$ abundance scale is on the right, the scale for the other species on the left.}
  \label{fig:crababun}
\end{figure}

This initial Crab model predicts that significant amounts of ArH$^+$ and OH$^+$ should be formed in gas of the observed composition and density subjected to the attenuated radiation field from the Crab PWN, in agreement with the observations of \citet{barlow2013}. This model predicts a higher abundance of OH$^+$ than ArH$^+$. However, \citet{barlow2013} found that the ArH$^+$ rotational line emission was significantly stronger than that from OH$^+$, so we calculated the predicted emissivities for both species. Figure~\ref{fig:crabemis} shows the line emissivities as a function of $A_V$ for the strongest transitions predicted by this Crab Nebula model, in the 450 to 1550 GHz range covered by the Herschel SPIRE spectrometer. The predicted emissivities of two [C I] and three OH$^+$ lines are higher than those of the two ArH$^+$ lines in this frequency range throughout the cloud, and the ArH$^+$ J=1-0 617 GHz transition is predicted to have a higher emissivity than the J=2-1 line at 1234 GHz. In contrast, \citet{barlow2013} found the J=2-1 line to be much stronger than the J=1-0 line, and of the other predicted lines only the OH$^+$ 971 GHz transition was detected, at a comparable strength to the ArH$^+$ 617 GHz line.

\begin{figure}
  \centering
  \includegraphics[width=\columnwidth]{fig4.eps}
  \caption{Emissivities of the strongest emission lines between 450 and 1550 GHz as a function of $A_V$ for the initial Crab Nebula model shown in Figure~\ref{fig:crababun}. ArH$^+$ transitions are shown as solid lines, OH$^+$ as dotted lines and neutral carbon as dashed lines.}
  \label{fig:crabemis}
\end{figure}

To investigate this discrepancy, a grid of models was run having hydrogen number densities of $n_{\rm{H}} = 1900$, $2 \times 10^4$, $2 \times 10^5$ and $2 \times 10^6$ cm$^{-3}$, and cosmic ray ionisation rates of $\zeta = 1$, $10^4$, $10^5$, $10^6$, $10^7$ and $10^8$ $\zeta_0$. These enhanced values cover the limits estimated by \citet{richardson2013}, based on the observed Crab synchrotron radiation and the degree of magnetic trapping. The grid was run using the UV and X-ray fluxes from the attenuated PWN spectrum described above. To compare surface brightnesses consistently between models with different densities, the line emissivities were integrated up to $A_V = 3$ in each model.

Figure~\ref{fig:linescrabd2e3} shows the surface brightness of the various transitions plotted in Figure~\ref{fig:crabemis}, relative to the OH$^+$ 971 GHz transition, with each ratio plotted as a function of $\zeta$, for models with $n_{\rm{H}} = 1900$ cm$^{-3}$. The dashed lines show the observed range of ratios relative to OH$^+$ 971 GHz for the two ArH$^+$ transitions. The relative surface brightnesses of the ArH$^+$ lines are much lower than observed, with the highest values occuring for $\zeta = 10^6 \zeta_0$. The ratio between the two ArH$^+$ transitions is also in conflict with observation except at the highest values of $\zeta$. The non-detection of the neutral carbon lines also favours a high value of $\zeta$, as carbon atoms are converted to C$^+$ ions by cosmic ray ionisation, while the relative intensities of the OH$^+$ lines are generally compatible with observations. Figure~\ref{fig:colinesd2e3} shows line surface brightnesses relative to the OH$^+$ 971 GHz transition versus $\zeta$ for the ten CO rotational transitions in the SPIRE FTS frequency range (J=4-3 to J=13-12). All the CO transitions are predicted to be significantly weaker than the OH$^+$ lines, in agreement with the lack of detected CO emission. With increasing $\zeta$, the rate of destruction of CO by cosmic rays and the temperature both rise, resulting in lower surface brightnesses (due to a lower CO abundance) and converging line strengths (as higher J levels can more easily be populated through collisions).

\begin{figure}
  \centering
  \includegraphics[width=\columnwidth]{fig5.eps}
  \caption{Line surface brightness ratios relative to OH$^+$ 971 GHz, of the other transitions of OH$^+$ (red), carbon (black) and ArH$^+$ J=1-0 (blue) and J=2-1 (green) from Figure~\ref{fig:crabemis} versus the cosmic ray ionisation rate $\zeta$, in units of the interstellar cosmic ray ionisation rate $\zeta_0 = 1.3 \times 10^{-17}$ s$^{-1}$, for models with $n_{\rm{H}} = 1900$ cm$^{-3}$. The dashed lines show upper and lower observational bounds for the ArH$^+$ transitions, from \citet{barlow2013}.}
  \label{fig:linescrabd2e3}
\end{figure}

\begin{figure}
  \centering
  \includegraphics[width=\columnwidth]{fig6.eps}
  \caption{CO line surface brightness ratios relative to OH$^+$ 971 GHz of rotational transitions from J=4-3 (upper) to J=13-12 (lower) versus the cosmic ray ionisation rate $\zeta$, in units of the interstellar cosmic ray ionisation rate $\zeta_0 = 1.3 \times 10^{-17}$ s$^{-1}$, for models with $n_{\rm{H}} = 1900$ cm$^{-3}$.}
  \label{fig:colinesd2e3}
\end{figure}

Figure~\ref{fig:linescrabd2e6} shows the line surface brightness ratios versus $\zeta$ for the same transitions as plotted in Figure~\ref{fig:linescrabd2e3}, for a hydrogen nucleus density of $2 \times 10^6$ cm$^{-3}$. None of the enhanced density models reproduce the observed ArH$^+$/OH$^+$ surface brightness ratios adequately, while for $n_{\rm{H}} = 2 \times 10^6$ cm$^{-3}$ even the highest $\zeta$ models still predict significant [C I] emission. The strength of the two ArH$^+$ lines relative to each other is best reproduced by models with high ($\ge 10^7 \zeta_0$) cosmic ray ionisation rates, which also produce the surface brightness ratios closest to those observed for ArH$^+$, and weak enough [C I] emission to be consistent with observation. All models predict OH$^+$ 1033 GHz emission approximately equal in strength to the SPIRE-FTS ArH$^+$ 617 GHz line. While this OH$^+$ transition is not widely detected, in one detector (SSW B1) we did find an emission line at 1033 GHz with a slightly lower surface brightness than that of the OH$^+$ 971 GHz line in the same detector, which would be consistent with these results.

\begin{figure}
  \centering
  \includegraphics[width=\columnwidth]{fig7.eps}
  \caption{Line surface brightness ratios relative to OH$^+$ 971 GHz versus the cosmic ray ionisation rate $\zeta$ in units of $\zeta_0$, for models with $n_{\rm{H}} = 2 \times 10^6$ cm$^{-3}$. The dashed lines show upper and lower observational bounds for the ArH$^+$ transitions, from \citet{barlow2013}.}
  \label{fig:linescrabd2e6}
\end{figure}

Figure~\ref{fig:colinesd2e6} shows the relative surface brightnesses of the CO transitions versus $\zeta$ for a cloud density of $2 \times 10^6$ cm$^{-3}$. The CO emission becomes stronger at higher densities, although it remains much weaker than OH$^+$ 971 GHz at high $\zeta$. The lack of detected CO emission excludes models with lower values of $\zeta$ and densities above $2 \times 10^4$ cm$^{-3}$, with the maximum value considered ($10^8 \zeta_0$) required to produce low enough CO emission at all densities.

\begin{figure}
  \centering
  \includegraphics[width=\columnwidth]{fig8.eps}
  \caption{CO line surface brightness ratios relative to OH$^+$ 971 GHz of rotational transitions from J=4-3 (upper) to J=13-12 (lower) versus the cosmic ray ionisation rate $\zeta$ in units of $\zeta_0$, for models with $n_{\rm{H}} = 2 \times 10^6$ cm$^{-3}$.}
  \label{fig:colinesd2e6}
\end{figure}

Figure~\ref{fig:oh971av3} shows the absolute surface brightness of the OH$^+$ 971 GHz transition as a function of $\zeta$, for the four densities investigated. \citet{barlow2013} found surface brightnesses of around $10^{-7}$ to $10^{-6}$ erg cm$^{-2}$ s$^{-1}$ sr$^{-1}$ for this line, much lower than the predicted values from our high $\zeta$ models, and for all our lowest density models. Figure~\ref{fig:oh971av1e-2} shows the surface brightness when the emissivity is integrated only up to a visual extinction of $A_V = 0.01$. In this case the $n_{\rm{H}} = 2 \times 10^4$ cm$^{-3}$ models have surface brightnesses closer to the observed values, with the three highest $\zeta$ models falling within the observed range. Figure~\ref{fig:linescrabd2e4av1e-2} shows line surface brightness ratios integrated up to $A_V = 0.01$ for a density of $2 \times 10^4$ cm$^{-3}$. Compared to the $A_V = 3$ case, the OH$^+$ and ArH$^+$ ratios are mostly unaffected, whereas there is a significant reduction in the predicted strength of the undetected [C I] transitions at low values of $\zeta$.

\begin{figure}
  \centering
  \includegraphics[width=\columnwidth]{fig9.eps}
  \caption{Absolute surface brightness integrated up to $A_V = 3$ of the OH$^+$ 971 GHz transition versus cosmic ray ionisation rate $\zeta$ for models with $n_{\rm{H}} = 1900$ cm$^{-3}$ (black), $2 \times 10^4$ cm$^{-3}$ (red), $2 \times 10^5$ cm$^{-3}$ (blue) and $2 \times 10^6$ cm$^{-3}$ (green). The dashed lines show upper and lower observational bounds from \citet{barlow2013}.}
  \label{fig:oh971av3}
\end{figure}

\begin{figure}
  \centering
  \includegraphics[width=\columnwidth]{fig10.eps}
  \caption{Absolute surface brightness integrated up to $A_V = 0.01$ of the OH$^+$ 971 GHz transition versus cosmic ray ionisation rate $\zeta$ for models with $n_{\rm{H}} = 1900$ cm$^{-3}$ (black), $2 \times 10^4$ cm$^{-3}$ (red), $2 \times 10^5$ cm$^{-3}$ (blue) and $2 \times 10^6$ cm$^{-3}$ (green). The dashed lines show upper and lower observational bounds from \citet{barlow2013}.}
  \label{fig:oh971av1e-2}
\end{figure}

\begin{figure}
  \centering
  \includegraphics[width=\columnwidth]{fig11.eps}
  \caption{Line surface brightness ratios relative to OH$^+$ 971 GHz versus the cosmic ray ionisation rate $\zeta$ in units of $\zeta_0$, for models with $n_{\rm{H}} = 2 \times 10^4$ cm$^{-3}$, with the emissivity integrated up to $A_V = 0.01$, as for our D4Z7 and D4Z8 models (see Table~\ref{tab:models}). The dashed lines show upper and lower observational bounds for the ArH$^+$ transitions, from \citet{barlow2013}.}
  \label{fig:linescrabd2e4av1e-2}
\end{figure}

A high value of $\zeta$ is required to account for the lack of neutral carbon and CO emission, but an increased X-ray/UV flux could also produce many of the same effects. To investigate this possibility, we repeated our \textsc{mocassin} modelling with the distance between the ionizing source and the filament reduced from $2.5$ pc to 1 pc. This resulted in input UV and X-ray field strengths of 94 Draines and 1.3 erg cm$^{-1}$ s$^{-1}$, respectively. Figure~\ref{fig:lines1pcd2e3} shows relative line surface brightnesses versus $\zeta$ for a cloud density of $1900$ cm$^{-3}$, using this increased radiation field. At low values of $\zeta$, the relative brightnesses of the neutral carbon lines are significantly reduced compared to the previous case, while the ArH$^+$ lines are increased slightly. However, to reduce the [C I] line surface brightnesses to levels consistent with observation, the cosmic ray ionization rate must still be at least $10^7 \zeta_0$, by which point the radiation field is no longer the dominant source of ionization and the two models give virtually identical results. We therefore conclude that an increased UV and X-ray radiation field cannot account for the observed emission line strengths, and high values of $\zeta$ are required.

\begin{figure}
  \centering
  \includegraphics[width=\columnwidth]{fig12.eps}
  \caption{Line surface brightness ratios relative to OH$^+$ 971 GHz versus the cosmic ray ionisation rate $\zeta$ in units of $\zeta_0$, for models with $n_{\rm{H}} = 1900$ cm$^{-3}$, for the Crab Nebula PWN spectrum at a distance of 1 pc. The dashed lines show upper and lower observational bounds for the ArH$^+$ transitions, from \citet{barlow2013}.}
  \label{fig:lines1pcd2e3}
\end{figure}

The rate coefficient used for the ArH$^+$ + $e^-$ reaction, $10^{-9}$ cm$^3$ s$^{-1}$, is an upper limit found experimentally by \citet{mitchell2005} for low energy collisions. As this upper limit is comparable to the rate coefficient for the ArH$^+$ + H$_2$ reaction ($8 \times 10^{-10}$ cm$^3$ s$^{-1}$), if the electron density is similar to that of molecular hydrogen the destruction rate of ArH$^+$ may be overestimated. The grid was rerun with reduced electron recombination rates of $10^{-10}$ and $10^{-11}$ cm$^3$ s$^{-1}$ in order to investigate the effects on the abundance of ArH$^+$. The highest gas temperature reached in any of the PDR models is around $12000$ K, corresponding to a thermal energy of about $1.0$ eV, while \citet{mitchell2005} found the reaction rate for electron energies below $2.5$ eV was negligible.

Figure~\ref{fig:lineserate-11d2e3} shows the relative line surface brightnesses versus $\zeta$ for a cloud density of $1900$ cm$^{-3}$, this time with the ArH$^+$ + e$^-$ reaction rate set to $10^{-11}$ cm$^3$ s$^{-1}$. For $\zeta \ge 10^6 \zeta_0$, the relative surface brightnesses of the two ArH$^+$ transitions have increased to values comparable with observations, although the absolute surface brightness of the OH$^+$ 971 GHz line is too high for these values of $\zeta$, as shown in Figure~\ref{fig:oh971av3}. Integrated only up to $A_V = 0.05$, the OH$^+$ surface brightness is within the observed range for all values of $\zeta$, while the relative strengths of the ArH$^+$ lines, shown in Figure~\ref{fig:lineserate-11d2e3av5e-2}, are broadly consistent with observation for $\zeta \ge 10^6 \zeta_0$, although for the $10^6 \zeta_0$ model the predicted [C I] 809 GHz surface brightness is too high. Figure~\ref{fig:lineserate-11d2e4av1e-2} shows the line surface brightness ratios integrated up to $A_V = 0.01$ for a density of $2 \times 10^4$ cm$^{-3}$, with the reduced dissociative recombination rate. Comparing with Figure~\ref{fig:linescrabd2e4av1e-2}, the ArH$^+$ line ratios are increased for $\zeta \ge 10^7 \zeta_0$. For these values of $\zeta$, all line surface brightnesses are consistent with observation, as both the relative line strengths and the OH$^+$ 971 GHz surface brightness fall within the observed range of values. The models with $n_{\rm{H}} = 2 \times 10^5$ and $2 \times 10^6$ cm$^{-3}$ do not reproduce the observed ArH$^+$ surface brightness ratios for any combination of $\zeta$ and $A_V$.

\begin{figure}
  \centering
  \includegraphics[width=\columnwidth]{fig13.eps}
  \caption{Line surface brightness ratios relative to OH$^+$ 971 GHz versus the cosmic ray ionisation rate $\zeta$ in units of $\zeta_0$, for models with $n_{\rm{H}} = 1900$ cm$^{-3}$ and a rate coefficient for ArH$^+$ + e$^-$ of $10^{-11}$ cm$^3$ s$^{-1}$. The dashed lines show upper and lower observational bounds for the ArH$^+$ transitions, from \citet{barlow2013}.}
  \label{fig:lineserate-11d2e3}
\end{figure}

\begin{figure}
  \centering
  \includegraphics[width=\columnwidth]{fig14.eps}
  \caption{Line surface brightness ratios relative to OH$^+$ 971 GHz versus the cosmic ray ionisation rate $\zeta$ in units of $\zeta_0$, for models with $n_{\rm{H}} = 1900$ cm$^{-3}$ and a rate coefficient for ArH$^+$ + e$^-$ of $10^{-11}$ cm$^3$ s$^{-1}$, with the emissivity integrated up to $A_V = 0.05$, as for our D3Z7 and D3Z8 models (see Table~\ref{tab:models}). The dashed lines show upper and lower observational bounds for the ArH$^+$ transitions, from \citet{barlow2013}.}
  \label{fig:lineserate-11d2e3av5e-2}
\end{figure}

\begin{figure}
  \centering
  \includegraphics[width=\columnwidth]{fig15.eps}
  \caption{Line surface brightness ratios relative to OH$^+$ 971 GHz versus the cosmic ray ionisation rate $\zeta$ in units of $\zeta_0$, for models with $n_{\rm{H}} = 2 \times 10^4$ cm$^{-3}$ and a rate coefficient for ArH$^+$ + e$^-$ of $10^{-11}$ cm$^3$ s$^{-1}$, with the emissivity integrated up to $A_V = 0.01$, as for our D4Z7 and D4Z8 models (see Table~\ref{tab:models}). The dashed lines show upper and lower observational bounds for the ArH$^+$ transitions, from \citet{barlow2013}.}
  \label{fig:lineserate-11d2e4av1e-2}
\end{figure}

\section{Discussion}
\label{sec:discussion}

Our initial models, and models with an increased UV/X-ray radiation field, were unable to reproduce the observed Herschel SPIRE FTS line emission for any combination of parameters we tested. With a reduced rate for the ArH$^+$ dissociative recombination reaction, results broadly consistent with observation can be found for densities of $1900$ and $2 \times 10^4$ cm$^{-3}$ and high cosmic ray ionization rates, with the highest density models ruled out. Cosmic ray ionization rates of at least $\zeta = 10^7 \zeta_0$ are required to reduce the [C I] emission to below detectable levels. Selecting filament thicknesses which reproduce the observed OH$^+$ 971 GHz surface brightnesses, we restrict our analysis to four models which are entirely consistent with the Herschel SPIRE FTS observations, D3Z7, D3Z8, D4Z7 and D4Z8, with the model parameters given in Table~\ref{tab:models}.

\subsection{Physical properties}

Figure~\ref{fig:properate-11d2e3z1e7} shows the gas temperature and abundances of H, H$_2$ and e$^-$ versus $A_V$ for an $n_{\rm{H}} = 1900$ cm$^{-3}$, $\zeta = 10^7 \zeta_0$ model, the values used for model D3Z7. Even for this case, with the lowest value of $\zeta$ of our preferred models, cosmic ray ionization dominates both the gas heating rate and the chemistry, so the filament properties do not change significantly with depth. The gas remains primarily atomic, with an ionized fraction of $\sim 0.4$ and a molecular fraction of $\sim 10^{-5}$, at a temperature of $\sim 8000$ K. For the $n_{\rm{H}} = 1900$ cm$^{-3}$, $\zeta = 10^8 \zeta_0$ model, there is no detectable change in any filament property with depth.

\begin{figure}
  \centering
  \includegraphics[width=\columnwidth]{fig16.eps}
  \caption{Temperature (red line) and abundance of H (solid black line), H$_2$ (dashed black line) and e$^-$ (dotted black line) relative to $n_{\rm{H}}$ versus $A_V$ for an $n_{\rm{H}} = 1900$ cm$^{-3}$, $\zeta = 10^7 \zeta_0$ model. These are the same parameters as model D3Z7, for which the cloud only extends up to $A_V = 0.05$.}
  \label{fig:properate-11d2e3z1e7}
\end{figure}

As the properties for our preferred models vary little with depth, we assume the values at the filament surface are representative for the whole filament and give the properties of our final models in Table~\ref{tab:prop}. We find that all models successfully reproducing the line surface brightness observations have H$_2$ abundances of order $10^{-5}$ and ionization fractions of at least $0.1$. We list the surface brightnesses of emission lines observed in the Crab Nebula for each model in Table~\ref{tab:emis}, along with observational upper and lower limits.

\begin{table*}
  \centering
  \caption{Gas temperatures and abundances of H, H$_2$, e$^-$, Ar$^+$ and ArH$^+$ for our final models.}
  \begin{tabular}{ccccccc}
    \hline
    Model & T / K & $n(\rm{H})$ / $n_{\rm{H}}$ & $n(\rm{H}_2)$ / $n_{\rm{H}}$ & $n(\rm{e}^-)$ / $n_{\rm{H}}$ & $n(\rm{Ar}^+)$ / $n_{\rm{H}}$ & $n(\rm{ArH}^+)$ / $n_{\rm{H}}$ \\
    \hline
    D3Z7 & $8381$ & $0.84$ & $2.1 \times 10^{-5}$ & $0.42$ & $6.5 \times 10^{-6}$ & $2.3 \times 10^{-8}$ \\
    D3Z8 & $11958$ & $0.53$ & $9.2 \times 10^{-6}$ & $1.24$ & $9.0 \times 10^{-6}$ & $8.1 \times 10^{-9}$ \\
    D4Z7 & $4745$ & $0.96$ & $5.0 \times 10^{-5}$ & $0.11$ & $2.2 \times 10^{-6}$ & $2.4 \times 10^{-8}$ \\
    D4Z8 & $8209$ & $0.84$ & $2.7 \times 10^{-5}$ & $0.42$ & $6.4 \times 10^{-6}$ & $2.9 \times 10^{-8}$ \\
    \hline
  \end{tabular}
  \label{tab:prop}
\end{table*}

\begin{table*}
  \centering
  \caption{Predicted surface brightnesses and observational limits for ArH$^+$, OH$^+$ and H$_2$ emission lines in the Crab Nebula, in erg cm$^{-2}$ s$^{-1}$ sr$^{-1}$. Observed values are taken from \citet{barlow2013} for ArH$^+$ and OH$^+$ and \citet{loh2011} for H$_2$.}
  \begin{tabular}{ccccc}
    \hline
    Model & OH$^+$ 971 GHz & ArH$^+$ J=1-0 & ArH$^+$ J=2-1 & H$_2$ $2.12$ $\mu$m \\
    \hline
    Observed & $3.4-10.3 \times 10^{-7}$ & $2.2-9.9 \times 10^{-7}$ & $1.0-3.8 \times 10^{-6}$ & $1.0-4.8 \times 10^{-5}$ \\
    D3Z7 & $6.0 \times 10^{-7}$ & $7.5 \times 10^{-7}$ & $9.9 \times 10^{-7}$ & $6.4 \times 10^{-7}$ \\
    D3Z8 & $6.2 \times 10^{-7}$ & $4.3 \times 10^{-7}$ & $1.0 \times 10^{-6}$ & $2.3 \times 10^{-7}$ \\
    D4Z7 & $4.1 \times 10^{-7}$ & $2.5 \times 10^{-7}$ & $5.6 \times 10^{-7}$ & $4.7 \times 10^{-7}$ \\
    D4Z8 & $7.8 \times 10^{-7}$ & $3.6 \times 10^{-7}$ & $2.1 \times 10^{-6}$ & $3.0 \times 10^{-7}$ \\
    AVD3Z7 & $1.0 \times 10^{-6}$ & $1.2 \times 10^{-6}$ & $1.7 \times 10^{-6}$ & $1.1 \times 10^{-6}$ \\
    AVD3Z8 & $1.0 \times 10^{-6}$ & $7.2 \times 10^{-7}$ & $1.7 \times 10^{-6}$ & $3.9 \times 10^{-7}$ \\
    AVD4Z7 & $3.4 \times 10^{-6}$ & $2.0 \times 10^{-6}$ & $4.7 \times 10^{-6}$ & $4.0 \times 10^{-6}$ \\
    AVD4Z8 & $6.5 \times 10^{-6}$ & $2.9 \times 10^{-6}$ & $1.8 \times 10^{-5}$ & $2.5 \times 10^{-6}$ \\
    \hline
  \end{tabular}
  \label{tab:emis}
\end{table*}

\subsection{H$_2$ vibration-rotation emission}

\citet{loh2010,loh2011} observed emission from the H$_2$ $2.12$ $\mu$m line from 55 knots in the Crab Nebula, with average surface brightnesses of $\sim 1-5 \times 10^{-5}$ erg cm$^{-2}$ s$^{-1}$ sr$^{-1}$. Figure~\ref{fig:h2212av3} shows the surface brightness of this line versus $\zeta$ for our models, integrated up to $A_V = 3$. Our $n_{\rm{H}} = 1900$ cm$^{-3}$ models give a 2.12 $\mu$m surface brightness consistent with observations, while for $n_{\rm{H}} = 2 \times 10^4$ cm$^{-3}$ the predicted surface brightnesses are somewhat higher than observed. However, these $A_V$ values give ArH$^+$ and OH$^+$ surface brightnesses far in excess of the observed values, as shown in Figure~\ref{fig:oh971av3}. For a filament thickness reduced to $A_V = 0.05$, shown in Figure~\ref{fig:h2212av5e-2}, the H$_2$ emission is reduced, with the models giving the highest values still an order of magnitude below observations. Our final models all predict H$_2$ surface brightnesses below the observed range, given in Table~\ref{tab:emis}.

\begin{figure}
  \centering
  \includegraphics[width=\columnwidth]{fig17.eps}
  \caption{Absolute surface brightness integrated up to $A_V = 3$ of the H$_2$ $2.12$ $\mu$m transition versus cosmic ray ionisation rate $\zeta$ for models with $n_{\rm{H}} = 1900$ cm$^{-3}$ (black), $2 \times 10^4$ cm$^{-3}$ (red), $2 \times 10^5$ cm$^{-3}$ (blue) and $2 \times 10^6$ cm$^{-3}$ (green). The dashed lines show upper and lower observational bounds from \citet{loh2011}.}
  \label{fig:h2212av3}
\end{figure}

\begin{figure}
  \centering
  \includegraphics[width=\columnwidth]{fig18.eps}
  \caption{Absolute surface brightness integrated up to $A_V = 0.05$ (as for models D3Z7 and D3Z8) of the H$_2$ $2.12$ $\mu$m transition versus cosmic ray ionisation rate $\zeta$ for models with $n_{\rm{H}} = 1900$ cm$^{-3}$ (black), $2 \times 10^4$ cm$^{-3}$ (red), $2 \times 10^5$ cm$^{-3}$ (blue) and $2 \times 10^6$ cm$^{-3}$ (green). The dashed lines show upper and lower observational bounds from \citet{loh2011}.}
  \label{fig:h2212av5e-2}
\end{figure}

\citet{loh2012} used H$_2$ vibration-rotation emission to derive the excitation temperature and density in Crab Nebula H$_2$ knots, finding a temperature range of $2000-3000$ K, and a lower limit of $n_{\rm{H}} = 2 \times 10^4$ cm$^{-3}$ for the density. While the limit on the density is consistent with two of our final models, D4Z7 and D4Z8, the gas temperatures for these models are higher than the derived excitation temperatures. Table~\ref{tab:h2lines} lists the surface brightnesses of the six H$_2$ vibrational transitions detected by \citet{loh2012}, relative to the 1-0 S(1) line, for a representative Crab Nebula knot from that paper and for our final PDR models. With the exception of D4Z8, all our models predict line ratios in agreement with the typical observed values. \citet{loh2012} also calculated the electron densities of the H$_2$ knots based on [S II] emission. Their typical electron abundances of $n(\rm{e}^-)/n_{\rm{H}} \approx 0.1$ are very close to the value of $0.11$ from our D4Z7 model, which also has the gas temperature closest to the measured H$_2$ exctiation temperature of all our final models.

\begin{table*}
  \centering
  \caption{H$_2$ vibrational line surface brightness relative to the 1-0 S(1) line, for Knot 51 from \citet{loh2012} and for our final models.}
  \begin{tabular}{ccccccc}
    \hline
    Line & 1-0 S(0) & 1-0 S(1) & 1-0 S(2) & 2-1 S(1) & 2-1 S(2) & 2-1 S(3) \\
    \hline
    Knot 51 & $0.23 \pm 0.04$ & $1.00 \pm 0.04$ & $0.52 \pm 0.09$ & $0.19 \pm 0.03$ & $<0.13$ & $<0.28$ \\
    D3Z7 & $0.16$ & $1.00$ & $0.44$ & $0.16$ & $0.06$ & $0.26$ \\
    D3Z8 & $0.16$ & $1.00$ & $0.44$ & $0.16$ & $0.06$ & $0.26$ \\
    D4Z7 & $0.18$ & $1.00$ & $0.43$ & $0.21$ & $0.09$ & $0.32$ \\
    D4Z8 & $0.18$ & $1.00$ & $0.48$ & $0.44$ & $0.20$ & $0.73$ \\
    \hline
  \end{tabular}
  \label{tab:h2lines}
\end{table*}

\citet{richardson2013} modelled a Crab Nebula H$_2$-emitting knot using \textsc{cloudy} \citep{ferland1998}, in order to explain the combination of near-infrared H$_2$ vibrational lines with optical emission from atomic and ionised gas. Their models required an additional heating source to reproduce the observed H$_2$ vibrational spectrum, which they assumed to be either energy input from shocks or turbulence, or from the ionising particle flux from the PWN. Their charged particle model corresponds to $\zeta = 6.1 \times 10^6 \zeta_0$, only slightly lower than the range of values we require. Their models assume the gas is in pressure balance, leading to a PDR region with the density increasing from $\sim 10^3$ to $10^5$ cm$^{-3}$. The gas temperature and H$_2$ and e$^-$ abundances are consistent with our final models where the densities overlap, although in the highest density regions the H$_2$ fractional abundance reaches $10^{-3}$, higher than any of our models. We note that in their models, these regions produce the majority of the H$_2$ $2.12$ $\mu$m emission, suggesting that an additional, denser gas component combined with our final models may solve the issue of underpredicting the H$_2$ surface brightness. However, as shown in Figure~\ref{fig:oh971av1e-2}, even a small amount of dense gas produces OH$^+$ 971 GHz emission in excess of the observed values, while Figure~\ref{fig:h2212av5e-2} shows that the H$_2$ $2.12$ $\mu$m surface brightnesses of these denser models are not significantly greater than for our $n_{\rm{H}} = 2 \times 10^4$ cm$^{-3}$ models.

\subsection{Far-infrared line emission}

\citet{gomez2012} listed Herschel PACS and ISO LWS fluxes for infrared fine structure emission lines from a number of species in the Crab Nebula, including [C II] and [O I], which are also modelled by \textsc{ucl\_pdr}. Table~\ref{tab:emis2} lists the predicted [O I] 63 $\mu$m/146 $\mu$m and [O I] 146 $\mu$m/[C II] 158 $\mu$m line ratios for our final models. The predicted ratios show reasonable agreement with the observed ratios from \citet{gomez2012}, although the D4Z7 and D4Z8 models predict [O I] 146 $\mu$m/[C II] 158 $\mu$m ratios twice as large as the highest value observed. However, we find that the predicted line surface brightnesses exceed the measured values by an order of magnitude. This overestimation may be due to the fact that our models do not treat any optical or UV line emission from higher ionization stages that would be produced by the high charged particle flux, and so require stronger emission from other species in order to balance the heating rate. [Fe II], which \citet{richardson2013} found to be the dominant coolant in the PDR region of their ionizing particles model, is also not included in our models.

\begin{table*}
  \centering
  \caption{Predicted and observed line surface brightness ratios for [O I], [C II] and HeH$^+$ emission lines in the Crab Nebula. Observed values are taken from \citet{gomez2012}.}
  \begin{tabular}{ccccc}
    \hline
    Model & $\frac{\rm{[O I]} 63 \mu\rm{m}}{\rm{[O I]} 146 \mu\rm{m}}$ & $\frac{\rm{[O I]} 146 \mu\rm{m}}{\rm{[C II]} 158 \mu\rm{m}}$ & $\frac{\rm{HeH}^+ 149 \mu\rm{m}}{\rm{[O I] 146} \mu\rm{m}}$ & $\frac{\rm{HeH}^+ 74.5 \mu\rm{m}}{\rm{[O I] 146} \mu\rm{m}}$ \\
    \hline
    Observed & $16.4-38.7$ & $0.125-0.323$ & - & - \\
    D3Z7 & $13.1$ & $0.303$ & $0.058$ & $0.055$ \\
    D3Z8 & $18.4$ & $0.213$ & $0.219$ & $0.221$ \\
    D4Z7 & $22.6$ & $0.650$ & $0.034$ & $0.034$ \\
    D4Z8 & $26.7$ & $0.582$ & $0.214$ & $0.250$ \\
    AVD3Z7 & $12.9$ & $0.305$ & $0.058$ & $0.055$ \\
    AVD3Z8 & $18.2$ & $0.214$ & $0.218$ & $0.220$ \\
    AVD4Z7 & $20.8$ & $0.660$ & $0.034$ & $0.034$ \\
    AVD4Z8 & $25.2$ & $0.585$ & $0.207$ & $0.247$ \\
    \hline
  \end{tabular}
  \label{tab:emis2}
\end{table*}

\subsection{Filament size and extinction}

The $A_V$ required to match the observed OH$^+$ 971 GHz surface brightness in our final models corresponds to filament thicknesses of $\sim 4 \times 10^{16}$ cm for $n_{\rm{H}} = 1900$ cm$^{-3}$ and $\sim 8 \times 10^{14}$ cm for $n_{\rm{H}} = 2 \times 10^4$ cm$^{-3}$. \citet{richardson2013} used a cloud thickness of $10^{16.5}$ cm, based on the observed width of the knot they modelled, close to the size of our $n_{\rm{H}} = 1900$ cm$^{-3}$ models. If the ionized region (which we do not attempt to model) is excluded, the PDR regions in their models are $\sim 5 \times 10^{15}$ cm in size, intermediate between our two values. \citet{loh2012} estimated the thickness of the H$_2$-emitting regions they observed to be of order $10^{14}$ cm, which is comparable to our $n_{\rm{H}} = 2 \times 10^4$ cm$^{-3}$ models.

\citet{grenman2017} used Hubble Space Telescope optical images of the Crab Nebula to determine the extinctions of several dusty globules observed as dark spots against the synchrotron background, obtaining values of $A_V = 0.20-0.34$. The globule sizes range from $6 \times 10^{15}$ cm to $3 \times 10^{16}$ cm, similar to the sizes of the H$_2$ emitting knots observed by \citet{loh2011}, although \citet{grenman2017} note that only $\sim 10 \%$ of their dusty globules are coincident with H$_2$ knots. While the observed globule sizes are similar to our $n_{\rm{H}} = 1900$ cm$^{-3}$ model cloud thicknesses, none of our preferred models match the extinctions found by \citet{grenman2017}. These models used a standard interstellar $A_V$/$N_{\rm{H}}$ ratio of $6.289 \times 10^{-22}$ cm$^2$ mag, whereas with the increased dust-to-gas ratio in the Crab Nebula the true value is likely to be higher.

We reran our final models with an $A_V$/$n_{\rm{H}}$ ratio increased by a factor of $3.6$ (the dust-to-gas mass ratio in the Crab relative to the standard ISM value) and an $A_V$ of $0.3$. The model parameters for these models, which we refer to as AVD3Z7, AVD3Z8, AVD4Z7 and AVD4Z8, are listed in Table~\ref{tab:models}. The temperature and abundances of the models were largely unchanged compared to models with the standard $A_V$/$n_{\rm{H}}$ and the same density and ionization rate, and the H$_2$ line ratios given in Table~\ref{tab:h2lines} were also essentially identical. We give the predicted line surface brightnesses in Tables~\ref{tab:emis} and \ref{tab:emis2}. The AVD4Z8 model predicts OH$^+$ and ArH$^+$ surface brightnesses significantly in excess of the observed values, while the other three models predict values either within or slightly above the observed range. The AVD4Z7 model additionally predicts an H$_2$ $2.12$ $\mu$m surface brightness only a factor of $2.5$ smaller than the observational lower limit, while also producing line ratios consistent with the observed values from \citet{loh2012}. The surface brightnesses of the [C II] and [O I] lines are higher than for the equivalent standard $A_V$/$n_{\rm{H}}$ models, while the line ratios are similar.

\subsection{Predicted HeH$^+$ rotational line emission}

HeH$^+$ has been predicted to form in ionized nebulae \citep{black1978,flower1979,roberge1982,cecchi1993}, as well as in the early Universe \citep{galli1998}, although the predicted intensities of both its rotational and vibrational transitions have been generally higher than the observational upper limts (e.g. \citet{moorhead1988,liu1997}). We incorporated HeH$^+$ into our chemical network using the reactions and rate coefficients from \citet{roberge1982}, with the exception of the He$^+$ + H reaction, for which we used the rate coefficient from \citet{zygelman1990}, and the HeH$^+$ + e$^-$ reaction, where we used the experimentally determined rate from \citet{yousif1989} (via \citet{cecchi1993}). We also neglected photodissociation of HeH$^+$, as our input SED from \textsc{mocassin} has negligible flux in the wavelength range beyond the Lyman limit relevant for the cross-section given by \citet{roberge1982}. We calculated the emissivity of the first three rotational transitions in the same way as for ArH$^+$, using electron impact excitation data from \citet{hamilton2016} and energy levels and A-coefficients from CDMS \citep{muller2001,muller2005}.

Our models predict HeH$^+$ J=1-0 $149.14$ $\mu$m and J=2-1 $74.57$ $\mu$m emission with similar surface brightnesses, from a few times $10^{-5}$ to a few times $10^{-4}$ erg cm$^{-2}$ s$^{-1}$ sr$^{-1}$, while the J=3-2 line is a factor of $\sim 20$ weaker. We give the line surface brightness ratios of the first two transitions compared to the [O I] 146 $\mu$m line in Table~\ref{tab:emis2} for our final models. For models with $\zeta = 10^7 \zeta_0$, the $149.14$ $\mu$m line is about 20 times weaker than the nearby [O I] 146 $\mu$m line, which is potentially detectable given the reported uncertainties in the [O I] line strengths, whereas for $10^8 \zeta_0$ the ratio increases to $\sim 0.2$ and the HeH$^+$ line should be clearly present. The $74.57$ $\mu$m ratios are similar to those for $149.14$ $\mu$m, or somewhat larger for the D4Z8 models, although the surface brightnesses are only around $1 \%$ of the value of the [O I] 63 $\mu$m line.

\subsection{Geometrical effects}

We treat the Crab Nebula filaments as homogenous one-dimensional slabs, with the PWN radiation field incident on one side and the observer on the other. In reality the Crab is a complex 3D structure, and all of these assumptions are likely to be incorrect to a greater or lesser extent. Each line of sight towards the Crab Nebula may intersect a number of gas components varying in density, temperature, radiation field strength and size, with the observed emission then being made up of multiple separate contributions. Additionally, the observed surface brightnesses are the average values over the area of the telescope beam, meaning that variations on smaller scales can be washed out. Our analysis, which assumes that the emission from all observed species originates from a single cloud, is obviously simplistic, but given the complexity of attempting to model multiple gas components using a 1D code, and our success at predicting molecular emission consistent with observation, we consider it acceptable for the purposes of this paper.

The effect of beam size on the observed surface brightnesses may resolve the issue of the excess [C II] and [O I] emission from our models. The dusty globules detected by \citet{grenman2017} have radii of $\leq 0.5$ arcsec, far smaller than the beam sizes in \citet{gomez2012} ($9.4\rm{''} \times 9.4\rm{''}$ for each PACS spaxel). The ratio of globule angular area to PACS spaxel area is $\sim 3 \times 10^{-3}$, which is low enough to explain the discrepancy between our models and the observations if each spaxel contains only a few globules, although unlike the moleuclar emission, there is no reason to assume that the far-infrared lines only originate in dense, shielded clumps. The Hershel SPIRE FTS beam size is also much larger than the typical angular sizes of the H$_2$ knots from \citet{loh2011}, so assuming that the ArH$^+$ and OH$^+$ emission also originates in these knots, a higher model surface brightness would be necessary to reproduce observed values. If this was obtained by increasing the filament $A_V$ value, it would also increase the H$_2$ $2.12$ $\mu$m surface brightness, potentially removing the issue of our models predicting lower values than those from \citet{loh2011}.

\subsection{Time-dependent effects}

Our PDR models assume that the gas is in equilibrium, with the chemistry and thermal balance both in steady state. However, as the Crab Nebula is a relatively young object ($\sim 950$ yr), this assumption is likely to be invalid. \citet{richardson2013} noted this issue and found that the timescale for the formation of H$^-$, which controls the formation of H$_2$, exceeds the age of the Crab in their models, indicating that time-dependent calculations are required. However, modelling the chemical evolution of the Crab Nebula from supernova to its current state is a complex problem well beyond the scope of this paper. We therefore focus on the timescales for reactions relevant to the formation of the molecules.

Our models find a minimum fractional ionization of Ar$^+$/Ar$^0$ $\approx 0.2$, while for lower densities and higher values of $\zeta$ the value increases. The main formation route for ArH$^+$, Ar$^+$ + H$_2$, has a rate of $\sim 10^{-9}$ cm$^3$ s$^{-1}$ \citep{rebrion1989}, so for an H$_2$ density of $\sim 1$ cm$^{-3}$ (for $n_{\rm{H}} = 2 \times 10^4$ cm$^{-3}$) the formation timescale is $\sim 10^9$ s, or around 30 years, much less than the age of the Crab. This suggests ArH$^+$ can form efficiently from Ar$^+$ under Crab Nebula conditions. However, this assumes an equilibrium abundance of H$_2$, which also has to form from the presumably initially fully ionized supernova ejecta. Using the values from our models, we find a similar timescale to \citet{richardson2013} of $\sim 10^4$ yr for the formation of H$_2$ by associative detachment, which is significantly longer than the age of the Crab. Although this suggests our value for the formation timescale of ArH$^+$ is an underestimate, we note that even with an order of magnitude increase our timescale is still shorter than the age of the Crab, and that earlier in the history of the Crab Nebula the gas was presumably both denser and more highly ionized, which would tend to shorten both formation timescales. Our models predict less H$_2$ emission than is actually observed in the Crab, raising the prospect that we may in fact be underestimating, rather than overestimating, the H$_2$ abundance. It therefore seems plausible that enough ArH$^+$ could have formed up to the present age of the Crab Nebula to account for the observed line emission.

In our chemical network, HeH$^+$ is formed by the reaction of He with H$^+$ and H$_2^+$, and the reaction of He$^+$ with H. The He$^+$/He ratio ranges from $0.03$ to $0.75$ depending on the model. Calculating the typical timescales of all these reactions (the time for He/He$^+$ to be converted to HeH$^+$), we find the dominant reaction is He$^+$ + H for all our preferred models, with a timescale shorter than the others by several orders of magnitude. However, even for this reaction the lowest formation timescale is at least a few $\times 10^3$ yr, and this value occurs for the model with the lowest He$^+$ abundance. It therefore appears that HeH$^+$ should form at a much slower rate and on a longer timescale than ArH$^+$.

\section{Conclusions}
\label{sec:conclusion}

We have modelled a Crab Nebula filament using a combination of photoionisation and PDR codes, varying the density and cosmic ray ionisation rate to create a grid of models. We find that a high energy particle ionisation rate of $>10^{7} \zeta_0$ is required in order to explain the lack of detected [C I] and CO emission in the Herschel SPIRE FTS spectra reported by \citet{barlow2013}, while the observed line surface brightness ratios of ArH$^+$ to OH$^+$ can only be reproduced with these high ionisation rates and a reduced electron recombination rate for ArH$^+$. This is consistent with both experimentally determined reaction rates and theoretical estimates of the charged particle flux in the Crab. The Herschel observations can be entirely explained by models with densities between $n_{\rm{H}} = 1900$ and $2 \times 10^4$ cm$^{-3}$, while models with higher densities fail to reproduce the observed ArH$^+$ to OH$^+$ line ratios. Our models agree with the conclusion of \citet{richardson2013} that the gas in the Crab Nebula filaments is primarily atomic, even in the regions with strong H$_2$ emission.

Using an increased $A_V$/$N_{\rm{H}}$ ratio more appropriate for the Crab Nebula, and a total extinction within the range found by \citet{grenman2017} for Crab Nebula dusty globules, we find that models with $n_{\rm{H}} = 1900$ cm$^{-3}$ reproduce the observed ArH$^+$ and OH$^+$ surface brightnesses measured by \citet{barlow2013}, while underpredicting the H$_2$ surface brightness compared to \citet{loh2010,loh2011,loh2012}. Our $n_{\rm{H}} = 2 \times 10^4$ cm$^{-3}$, $\zeta = 10^7 \zeta_0$ model produces surface brightnesses within a factor of a few of observations for ArH$^+$, OH$^+$ and H$_2$, while the $n_{\rm{H}} = 2 \times 10^4$ cm$^{-3}$, $\zeta = 10^8 \zeta_0$ model predicts significantly stronger ArH$^+$ and OH$^+$ emission, as well as failing to reproduce the H$_2$ vibration-rotation line ratios. All models predict [C II] and [O I] line emission in excess of the values listed by \citet{gomez2012}, although surface brightness ratios are comparable with observed values. We suggest that these dusty globules are the source of most of the observed molecular emission in the Crab Nebula, with denser gas accounting for the $\sim 10 \%$ of globules associated with H$_2$ vibrational emission. Although our models predict HeH$^+$ rotational emission above detection thresholds, the formation timescales for this molecule are much longer than the age of the Crab Nebula, lowering the expected HeH$^+$ line strengths. The formation timescale for ArH$^+$, in contrast, is significantly shorter, suggesting that it can be formed efficiently under the conditions prevalent in the Crab Nebula.

\section*{Acknowledgements}

FP is supported by the Perren fund and IMPACT fund. The authors would like to acknowledge the work of Dr. Tom Bell, who was responsible for rewriting and substantially upgrading \textsc{ucl\_pdr}.

\bibliographystyle{mnras}