\section{Introduction}
Coreference resolution is the task of identifying all expressions in text that refer to the same entity.  In this paper we set out to provide an in-depth analysis and improve performance on a special coreference subtask: finding all references---either by name, pronoun or nominal---to a person named in the text. Named Person Coreference (NPC) would be especially useful for downstream information extraction tasks, in which facts about the person are extracted from large textual corpora and used as knowledge source for a range of artificial intelligence systems.  NPC extends the power of Named Entity Recognition (NER) systems by not only finding all text snippets that are a person's name, but also identifying all places in the text where the same person was mentioned.

Our work is oriented towards practical uses of the results in downstream applications, so we re-examine standard coreference metrics and find them lacking in their ability to quantify the performance of systems. Since applications require information about people and people are identified by their names, the evaluation metrics for this NPC task should focus on the resolution of mentions to the correct name. If all the pronouns referring to a person are resolved correctly to each other but are not linked to any named mention or are linked to a wrong named mention, their correct resolution to each other would not be useful for downstream applications. Standard coreference metrics do not incorporate these aspects of performance and hence give high performance for results unsuitable for further use. We also show that the existing metrics are not sensitive to finding any mention to a person at all. They give higher values for systems that do not find a large number of entities but do good coreference resolution on the subset of entities they find.

We introduce new metrics to overcome these shortcomings. We separate the identification of entities and resolution of different mention types, thus transparently tracking areas of system performance and improvement. We create a subset of the standard coreference resolution annotated data sets by identifying all named person entities and adapting existing coreference systems to filter their original output to only person entities identified by name.

Inspired by the error analysis of the performance of current coreference systems, we introduce a new solution to the NPC task that builds upon the capabilities of state of the art named entity recognition systems. We describe a highly effective rule-based approach to named person coreference that combines NER-driven rules and basic heuristics based on writing styles for resolving pronominal coreference. We show that this system has comparable performance to the state of the art systems on the standard metrics and superior performance on the task-informed metrics we introduce.

Finally, we test our approach for named entity coreference on entities other than people. We find that the state-of-the-art coreference systems do not suffer from the same issues for other entities as they do for persons and thus the results are similar across both standard metrics and NPC metrics.

\section{Why Named Person Coreference?}

People, and the information about people expressed in text, are of interest in numerous applications. Moreover, references to people behave differently from references to other named entities, indicating that the intersection of coreference and mentions to people represents a fruitful area for application-oriented research. We expand on these observations motivating our work.

\subsection{NPC in Downstream Applications}

Many information extraction and language technology tasks involve people. About 15\% of all web searches contain a person's name \cite{DBLP:conf/sigir/WeerkampBKMBR11} and especially in news search, NPC can help find articles in which the person of interest is the focus of discussion, mentioned by pronouns in addition to their full name.\footnote{An open problem in search is to disambiguate between different people with the same name, mentioned in different documents \cite{Artiles2010WePS3EC}. We do not deal with this problem and resolve mentions only within the same document.}

People are also often targets for information extraction systems \cite{Ji:2011:KBP:2002472.2002618} and for knowledge base completion tasks \cite{west2014knowledge}. Yet about half of the references to a person in text are not by name (cf. the first column of Table \ref{StatsOtherEntities}). Systems need to extract information about the person where they are mentioned explicitly by name as well as referenced by a pronoun.

Biography summarization \cite{zhou2005multi} needs to extract sentences from the news that contain information about a given person. More relevant sentences can be extracted if we know which pronouns and nominals refer to this person. Similarly, creation of proper noun ontologies \cite{mann2002fine} can use patterns other than (proper noun - common noun) if other references to the entity are known.

\ignore{There has also been interest in coreference resolution for person entities specifically in the context of video transcripts \cite{ramanathan2014linking}. }

\begin{table}[t]
\centering
\small
\setlength{\tabcolsep}{4pt}
\begin{tabularx}{\linewidth}{L{1.5}C{0.4}C{0.4}C{0.4}C{0.4}}
\hline
&{\bf PER}&{\bf ORG}&{\bf GPE}&{\bf DATE}\\\hline
Non-singleton &67.98\%&50.76\%&52.02\%&21.3\%\\
Entities &15.44\%&11.35\%&11.36\%&4.24\%\\
Entity Mentions&24.47\%&14.33\%&13.62\%&3.41\%\\
Named mentions&45.12\%&55.79\%&68.43\%&82.69\%\\
Non-named mentions&54.88\%&44.21\%&31.57\%&17.31\%\\
Avg cluster size&5.74&4.57&4.34&2.89\\\hline
\end{tabularx}
\caption{Statistics, in OntoNotes (nw,bn,mz), on coreference for PERson, ORGanization, geopolitical entity (GPE) and DATE named entities.  Non-singleton entities are mentioned at least twice in a text and so require coreference. Entities is the percentage of coreference clusters to entities of the given type. Entity mentions is the percentage of all individual references to any entity of a given type. Average cluster size is the average number of coreferent mentions.}\label{StatsOtherEntities}
\end{table}

\subsection{Coreference and Named Entities}
References to people are distributed quite differently from other types of named entities, as we show in Table \ref{StatsOtherEntities}. To compile these numbers, we make use of the OntoNotes coreference resolution corpus \cite{pradhan2007ontonotes} and gold-standard annotations for named entity recognition on the same data. In this way, we can quantify the patterns in coreference of different named entity types.

People are usually identified by their name in text:  86\% of the coreference clusters on people in the OntoNotes training data have at least one named mention; 88.2\% of animate third person singular pronouns (pronouns that can be used to refer to people) are a part of coreference chains with a named person, as seen in Table ~\ref{PronounStats}.\footnote{Female pronouns occur considerably less often than male pronouns, and a larger portion of female pronouns are to unnamed women. These numbers raise question about usage but overall the trend is clear: third person animate pronouns, in their vast majority, resolve to a person named in the article.} Therefore, it is uncommon to refer to a person without identifying them by their name at least once.

All named entities are on average much less likely to be singletons than a typical entity, mentioned only once in the text and not requiring coreference resolution \cite{DeMarneffe:2015:MLD:2831407.2831417}, and of the named entities people are most likely to be mentioned repeatedly. People are rarely mentioned only once in news (see the first line in Table \ref{StatsOtherEntities}): 68\% of people named in text have at least one other coreferent mention to them in contrast to 51\% and 52\% for organizations and locations respectively.

Named people entities make up 15\% of all coreference clusters in OntoNotes (see line 2 of Table \ref{StatsOtherEntities}), yet 25\% of all mentions that require coreference resolution are mentions of people (Line 3 in Table \ref{StatsOtherEntities}). There are more mentions of the same person on average (line 6 in Table \ref{StatsOtherEntities}) and PERson is the entity type with the largest portion of references that are not by name (line 5 in Table \ref{StatsOtherEntities}). Only 45\% of mentions to people are by name, compared to 56\% for organizations and 68\% for geo-political entities and over 80\% for dates.

In sum, references to named people make up a quarter of all references involved in coreference resolution and only half of these references are by name, so output of traditional NER systems that identify the names are not sufficient to track all mentions of the person; rules to disambiguate pronouns and noun phrase references are needed to track all mentions of people.

\begin{table}[t]
\centering
\small
\setlength{\tabcolsep}{4pt}
\begin{tabularx}{0.7\linewidth}{C{0.25}|C{0.25}|C{0.4}}
&{\bf Total}&{\bf Part of PER cluster}\\\hline
he&3062&2705 (88.3\%)\\
him&447&379 (84.7\%)\\
his&1931&1669 (86.4\%)\\\hline
she&718&588 (81.8\%)\\
her&667&521 (78.1\%)\\
hers&2&2 (100\%)\\\hline
\end{tabularx}
\caption{Third person singular pronouns in the Ontonotes train set.}\label{PronounStats}
\end{table}

\section{NPC and Coreference Evaluation}
\label{sec:generalmetrics}

Named Person Coreference is motivated by the needs in downstream applications. This setting allows us to critically review current practices for coreference evaluation and identify aspects where they fall short in quantifying system performance.

Shared tasks on coreference (at CoNLL-2011 and 2012 \cite{pradhan-EtAl:2014:P14-2} ) use the average of three scores as their official evaluation: MUC \cite{Vilain:1995:MCS:1072399.1072405}, \(B^3\) \cite{DBLP:conf/acl/BaggaB98} and CEAFE \cite{DBLP:conf/naacl/Luo05}.  \ignore{The first two scores,  MUC \cite{Vilain:1995:MCS:1072399.1072405} and \(B^3\) \cite{DBLP:conf/acl/BaggaB98}, assess system performance in terms of their ability to decide if a pair of mentions (noun phrases) in text refer to the same entity and  mention-focused recall/precision respectively. Below we describe these in more detail and point out how they are deficient when examined in the context of downstream tasks, particularly ones that involve NPC. The third score, CEAFE \cite{DBLP:conf/naacl/Luo05}, instead incorporates the notion of an entity as the key element in the evaluation.} Prior work \cite{moosavi2016coreference} discussed the shortcoming of these metrics and introduced the link entity aware (LEA) score to tackle them. Below we describe these and point out how they are deficient when examined in the context of downstream tasks, particularly ones that involve NPC. \footnote{This view of evaluation focused on the ability of a system to discover and track entities, which we later expand on, is similar to the vastly popular entity linking task, in which named mentions in text are linked to an abstract entity, such as one defined in Wikipedia \cite{Mihalcea:2007:WLD:1321440.1321475,DBLP:conf/acl/HanS11,DBLP:journals/tacl/DurrettK14,DBLP:conf/naacl/PanCHJK15,Radhakrishnan2018ELDENIE}.}

\begin{table}[t]
\centering
\small
\setlength{\tabcolsep}{4pt}
\begin{tabularx}{\linewidth}{L{}}
\\\hline
Gold clusters: \(\{John Doe, he_5, he_6, he_7, he_8, he_9\},\)\\
\(\{Richard Roe, he_1, he_2\}, \{Joe Smith, he_3, he_4\},\)\\\hline
Solution 1: \(\{John Doe, he_5, he_6, he_7, he_8, he_9\}\)\\
Solution 2: \(\{he_1, he_2\}, \{he_3, he_4\}, \{he_5, he_6, he_7, he_8, he_9\}\)
Solution 3: \(\{Richard Roe, he_1\}, \{Joe Smith, he_3\},\)\\
\(\{John Doe, he_5, he_6\}\)\\\hline
\end{tabularx}
\caption{NPC Examples}\label{Examples}
\end{table}

\paragraph{MUC}
MUC computes performance at the entity level, where a goldstandard cluster of mentions (noun phrases in text) represents an entity. The recall for an entity is based on
the minimum number of links that would have to be added in the predicted clusters containing any mention of this entity, to make them connected and part of the same cluster. Precision is computed by reversing the role of gold and predicted clusters.

\paragraph{B-cubed}
\(B^3\) measures performance on the mention level. It iterates over all goldstandard mentions of an entity, averaging the recall of its gold cluster in its predicted cluster. Different mentions may have been predicted as referring to different entities, placed in different solution clusters. It computes precision by reversing the role of gold and predicted clusters.

\paragraph{CEAF}
CEAF first finds a one-to-one mapping between gold and predicted clusters, in effect placing the discovery of entities first. It then computes recall as the number of same or similar mentions shared by the gold and predicted clusters divided by the number of mentions in the gold cluster. Precision similarly is equal to the number of same/similar entities in the gold and predicted, divided by the number of mentions in the predicted cluster. The final version averages these number either per mention as in B-cubed (CEAFm), or per entity (CEAFe), as would be more reasonable.

\paragraph{LEA}
LEA is a link-based metric most similar to MUC, that computes recall as the number of correctly resolved links between mentions, weighting the results for each entity by its number of mentions. In this way, resolving correctly an entity with more mentions contributes more to the overall score than resolving correctly all mentions to an entity mentions only twice for example. Precision is computed by reversing the role of gold and predicted clusters.

\begin{table}[t]
\centering
\small
\setlength{\tabcolsep}{2pt}
\begin{tabularx}{\linewidth}{C{0.35}|C{0.13}C{0.13}C{0.13}|C{0.13}C{0.13}C{0.13}|C{0.13}C{0.13}C{0.13}}
&\multicolumn{3}{>{\hsize=\dimexpr3\hsize+4\tabcolsep+1\arrayrulewidth\relax}C{0.13}|}{\bf Solution 1}&\multicolumn{3}{>{\hsize=\dimexpr3\hsize+4\tabcolsep+1\arrayrulewidth\relax}C{0.13}|}{\bf Solution 2}&\multicolumn{3}{>{\hsize=\dimexpr3\hsize+4\tabcolsep+1\arrayrulewidth\relax}C{0.13}}{\bf Solution 3}\\\cline{2-10}
&{\bf R}&{\bf P}&{\bf F1}&{\bf R}&{\bf P}&{\bf F1}&{\bf R}&{\bf P}&{\bf F1}\\\hline
MUC&0.55&1&0.71&0.66&1&0.8&0.44&1&0.61\\
B-cub&0.5&1&0.66&0.56&1&0.72&0.34&1&0.51\\
CEAFm&0.5&1&0.66&0.75&1&0.85&0.58&1&0.73\\
CEAFe&0.33&1&0.5&0.83&0.83&0.83&0.75&0.75&0.75\\
LEA&0.5&1&0.66&0.5&1&0.66&0.26&1&0.42
\end{tabularx}
\caption{Evaluation of the hypothetical solutions on NPC examples in Table \ref{Examples}}\label{ExampleStats}
\end{table}

Recall that the goal of NPC is to find all mentions referring to a person, identified by their name. We provide a made up example of a goldstandard and three possible solutions in Table ~\ref{Examples}. The gold standard contains three entities: John Doe, Richard Roe and Joe Smith. As the clusters indicate, each is also mentioned by a pronoun a number of times. The first hypothetical system response (solution 1) identifies only one entity.\ignore{, decides that there is only one entity mentioned in the text.} It finds all mentions to John Doe correctly but completely misses all mentions to the other two entities. The second solution correctly resolves the pronouns to each other but does not link these to any names. The third solution is able to identify a few mentions for all of the three entities.
Intuitively, solution 2 has little practical value, solution 3 is best and solution 1 is acceptable.

The results for standard coreference scores are shown in Table ~\ref{ExampleStats}.  All of the metrics have the highest values for the Solution 2, which does not identify a single name. The same would be true even if any of these correct set of pronouns were linked to the wrong name. This is because none of the metrics take into account the types of mentions and the need for resolution to correct names.

\section{NPC evaluation metrics}
\label{sec:npcmetrics}
We showed in the previous section that the existing coreference metrics are not suitable for the NPC task. Moreover, having a single number to track all aspects of the systems, makes it less interpretable. Prior work has argued that when coreference is used in downstream applications, evaluation criteria should be interpretable \cite{tuggener2014coreference}. We concur, and introduce a set of task-specific criteria for evaluation of NPC. These are inspired by error analysis we performed. For clarity, we provide examples errors on Ontonotes.

\subsection{Entity F1: Matching Output to Goldstandard Entities}

In the goldstandard, all noun phrases referring to the same person mentioned by name at least once are grouped into entities. In the system output, we also wish to find the chain corresponding to each person, similar to the motivation of the CEAFE evaluation. To map entity chains between the goldstandard and the system output, we select for each goldstandard entity, the predicted entity that has the highest F1 score with respect to the mentions it contains. The entity F1 score is the harmonic mean of the precision\footnote{Number of mentions to the goldstandard entity also in the system entity divided by the number of all system mentions.} and recall\footnote{Number of mentions to the goldstandard entity also in the system entity divided by the number of mentions in the goldstandard entity.} for the entity mapping.

To compute the intersection between a goldstandard and a system entity, we first augment each goldstandard entity with a list of all variations of the person's name. We rely on the goldstandard named entity annotation in OntoNotes and intersect this with the membership in a coreference chain. This provides lists of the full name, last name, occasionally nicknames and variants of the name, i.e. {\em \{Frank Curzio, Francis X. Curzio, Curzio\}}, {\em \{Dwayne Dog Chapman, Dog Chapman, Chapman\}}. We consider a predicted chain to be a candidate match for a gold chain only if it contains at least one of the name variants.

We consider a predicted mention to match if it matches the gold mention exactly in the same place of the text. Mention detection has a huge impact on evaluation but we use an exact span matching to be consistent with the existing coreference scorers.\footnote{Exact mention match is used for calculating F1. We do not use exact mention match to find candidate chains as the presence of the name can indicate which person the cluster is about.}

If a goldstandard entity does not get paired with any system entity, the F1 for that entity is taken to be zero. In our example above, Solution 1 will have poor recall, because it finds only one entity in a document containing three entities. We find the overall F1 of the system as the average of the F1 for each gold entity.

\subsection{Entity not found}

The entity F1 evaluation gives a sense of overall system performance but mixes true purity of the system-discovered entities and the ability to discover entities at all.
"Entity not found" is the error when no NPC system output overlaps with a gold standard entity. Entities not found contribute a score of 0 for the average F1.\footnote{We consider only chains containing a named mention. Chains that do not contain any named mention are filtered out. More details on filtering to follow in the section on system performance.}

\subsection{Pronoun Resolution Accuracy}
Finally, we track the entity F1 when only mentions of given syntactic type are preserved in the chain---name, pronoun and nominal. Of special interest is to track system performance when resolving pronouns. Many of these issues arise due to the need for commonsense knowledge and reasoning for correct resolution, as in the these examples:

\noindent
\textit{\(\textbf{A school official}_1\) talked about the Vice President's chances during an interview with the Boston Globe. \(\textbf{He}_2\) says it's unlikely Gore will be selected, because \(\textbf{he}_3\) doesn't have enough experience in the academic world.}

The pronoun \(\textit{he}_3\) is incorrectly resolved to the school official instead of Al Gore. To correctly resolve this, we need to know that not having enough experience would be the reason for Gore not getting selected and not the reason for the school official making a statement.

\noindent
\textit{Maybe Lily became so obsessed with where people slept and how because her own arrangements kept shifting. When \(\textbf{Rosie}_1\) died, \(\textbf{her}_2\) uncles moved in and let \(\textbf{her}_3\) make the sleeping and other household arrangements.}

Rosie is Lily's mother (explained earlier in the article). Both \(\textit{her}_{2,3}\) are incorrectly resolved to Rosie, despite the common sense fact that a person can not make household arrangements after their death.

\begin{table*}[t]
\centering
\small
\setlength{\tabcolsep}{4pt}
\begin{tabularx}{\linewidth}{L{1.2}C{0.31}C{0.31}C{0.31}C{0.31}|C{0.31}C{0.35}C{0.35}|C{0.28}C{0.28}C{0.28}}
\hline
&{\bf Chains not found}&{\bf F1 (NPC)}&{\bf Avg F1 (coref)}&{\bf LEA F1}&{\bf F1 (names)}&{\bf F1 (pronouns)}&{\bf F1 (nominals)}&{\bf F1 (nw)}&{\bf F1 (bn)}&{\bf F1 (mz)}\\\hline
CoreNlp deterministic&{\bf 24.76\%}&0.501&0.499&0.409&0.489&0.384&0.027&0.499&0.494&0.525\\
CoreNlp statistical&39.30\%&0.45&0.577&0.499&0.488&0.307&0.019&0.411&0.49&0.45\\
CoreNlp neural&27.35\%&{\bf 0.572}&0.677&0.613&{\bf 0.612}&{\bf 0.397}&0.062&{\bf 0.59}&0.58&0.502\\
AllenNlp&31.87\%&0.563&{\bf 0.71}&{\bf 0.66}&0.61&0.349&{\bf 0.078}&0.412&{\bf 0.699}&{\bf 0.608}\\\hline
\end{tabularx}
\caption{Performance of existing systems. The left panel shows named people coreference metrics (percentage chains not found and entity F1), and average F1 coreference evaluation combining MUC, \(B^3\) and CEAFE, on all test data. The middle panel shows entity F1 by type of mention, names, pronouns or nominals. The right panel shows entity F1 broken down by subgenres in the test data: newswire (nw), broadcast news (bn) and magazines (mz). Existing F1 metric is not sensitive to chains not found, and this explains the difference from the new intuitive NPC metric. For each metric, the best system has been bold-faced}\label{AllPeopleMetrics}
\end{table*}

\subsection{Over-splitting and over-combination of entities}
Sometimes, systems produce more than one clusters, each containing the same name. In this case we say that the gold-standard entity was over-split. Following is an example illustrating this error -

\noindent
\textit{\(\textbf{Frank Curzio}_1\). Many people now claim to have predicted the 1987 crash. \(\textbf{Queens newsletter writer Francis X. Curzio}_2\) actually did it: \(\textbf{He}_3\) stated in writing in September 1987 that the Dow Jones Industrial Average was likely to decline about 500 points the following month. \(\textbf{Mr. Curzio}_4\) says what happens now will depend a good deal on the Federal Reserve Board. If it promptly cuts the discount rate it charges on loans to banks, \(\textbf{he}_5\) says, '' That could quiet things down. ''}

The above text mentions Frank Curzio 5 times, however all the mentions aren't considered to be coreferent to each other. Mentions 2 and 3 form one cluster while the other three form another one.

At times systems also produce coreference chains that combines mentions to two different people. This error occurs when different people mentioned in the text have the same last name (notoriously for U.S. news Bill Clinton and Hilary Clinton) but also occasionally in cases where the names are completely different but the roles of the people are similar, as in the example below -

\noindent
\textit{\textbf{UN Secretary General Kofi \(\textbf {Annan}_1\)} said Wednesday, it is important to help spread democracy around the world.
VOA's Breck Ardery reports from the United Nations . In a new report, \(\textbf{the Secretary General}_2\) says democratization has now taken root as a universal norm and that the United Nations should strengthen its commitment to assisting nations that are moving toward democracy. Commenting on the report, \textbf{UN Assistant Secretary General for Political Affairs Danilo \(\textbf {Turk}_3\)} told reporters the principle of national sovereignty does not preclude support for democracy.}

Although it is clear from the above text snippet that Kofi Annan and Danilo Turk are two different people, they have similar titles and recognized to be coreferent.

Following is another document, where Scott Peterson and his wife Laci Peterson are said to be coreferrent, in spite of having distinct gender pronouns referring to each. \ignore{Apart from the highlighted mention, Scott Peterson in the second to last sentence isn't detected, and 'Peterson' and 'his wife Laci' in the last sentence form singleton clusters each.}

\noindent
\textit{Moving on to the \(\textbf{Laci Peterson}_1\) case Bill, an unusual request from \(\textbf{Scott Peterson}_2\)'s attorney. A pair of \textbf{Laci \(\textbf {Peterson}_3\)}'s missing shoes could be very important evidence in \(\textbf{Scott's}_4\) murder trial. \(\textbf{She}_5\) is asking anyone who finds the pair to give them back. No sign of Geragos or \(\textbf{Peterson}_6\) in court yesterday when a judge was considering whether to unseal warrants obtained before Scott Peterson's arrest. Peterson awaiting trial in the murder of \(\textbf{his}_7\) wife Laci and their unborn son.}

\section{Evaluation of existing systems}
\label{sec:evals}

We evaluate the Stanford coreference system, with its deterministic \cite{DBLP:conf/emnlp/RaghunathanLRCSJM10}, statistical \cite{DBLP:conf/acl/ClarkM15} and neural \cite{DBLP:conf/emnlp/ClarkM16} versions, and the neural end-to-end AllenNLP system \cite{DBLP:conf/emnlp/LeeHLZ17} on both the existing metrics and the NPC metrics.

These general coreference systems find coreferring expressions of any type and produce coreference clusters for all mentioned entities (groups  of noun phrases that refer to the same entity). In named person coreference, the goal is to find all mentions to {\em a person} who has been {\em referred to by name} at least once in the document. This means that the output of off-the-shelf coreference systems has to be filtered to keep only chains that contain at least one mention noun phrase with a syntactic head that is a person's name. For our evaluation, we use dependency parsing to detect whether a name is the head of a mention, by checking that no other word in the mention is an ancestor of the name in the dependency parse tree. We use automatic {\sc person} NER tags from Stanford coreNLP  \cite{DBLP:conf/acl/FinkelGM05} to determine if the head is a name. We call the coreference chains remaining after filtering {\em NPC entities}.

Less strict filtering, such as the presence of 3rd person singular personal or possessive pronouns would also indicate that the corresponding entity is a person. For NPC, we insist on having at least one named mention,\ignore{to facilitate search and information extraction,} as elaborated upon in the earlier sections.  We found that the AllenNLP system does not have a named mention in about 30\% of the coreference chains that do contain a personal or possessive third person pronoun. This number is about 20\% for the CoreNLP neural system.

For evaluating the systems, we use only the relevant subsample of the standard coreference evaluation data in OntoNotes. We work with the test newswire (89 documents, 136 NPC chains), broadcast news (93 documents, 135 NPC chains) and magazine (45 documents, 49 NPC chains) documents only. We evaluate on the 320 NPC chains in these 227 documents. The NPC chains contain three types of mentions---names, pronouns, and nominals. Nominals account for less than 5\% of the mentions in all genres, while the remaining mentions are split almost equally between  names and pronouns.

The third column in the first panel of Table~\ref{AllPeopleMetrics} shows the standard F1 on all the systems.  As expected, the AllenNLP systems outperforms all the systems, with CoreNLP neural as a close second. Both perform much than the CoreNLP deterministic and statistical systems.

However, this difference in performance isn't as big on the NPC F1 (see second column of Table \ref{AllPeopleMetrics}). The Stanford CoreNLP neural system has the highest NPC F1, only slightly better than the AllenNLP system.\footnote{This evaluation uses an exact span matching to be consistent with the existing coreference scorers. Our experiments showed a relaxed mention span matching allowing a difference of a few words results in gains up to 12 points F1. The gain with relaxed mention matching is highest for the deterministic system which ends up performing as good as the CoreNLP neural system.} The NPC F1 includes 0s for goldstandard entities not found by the system. We track this separately as well (first column of Table \ref{AllPeopleMetrics}). The deterministic and neural Stanford CoreNLP systems have the lowest percentage of entities not found errors, about a quarter of all entities. The Stanford statistical system is worst at finding entities, missing almost 40\% of entities.

We also separate the performance of the systems by mention type. The second panel of Table~\ref{AllPeopleMetrics} reveals that the systems make more mistakes on pronouns compared to names.

We also evaluated separately on the three genres, as shown in the third panel of Table \ref{AllPeopleMetrics}. While the CoreNLP system outperforms on newswire by a huge margin, the AllenNLP system performs the best on broadcast news and magazines. Like overall F1, scores are greatly impacted by the mention detection and the difference in performance becomes less with a relaxed mention matching.

We tracked the over-splitting and the over-combination of entities as well. However, the overall values were quite small and similar for all the systems and have thus not been included in the results here.

\begin{table*}[t]
\centering
\small
\setlength{\tabcolsep}{4pt}
\begin{tabularx}{\linewidth}{L{1.2}C{0.31}C{0.31}C{0.31}C{0.31}|C{0.31}C{0.35}C{0.35}|C{0.28}C{0.28}C{0.28}}
\hline
&{\bf Chains not found}&{\bf F1 (NPC)}&{\bf Avg F1 (coref)}&{\bf LEA F1}&{\bf F1 (names)}&{\bf F1 (pronouns)}&{\bf F1 (nominals)}&{\bf F1 (nw)}&{\bf F1 (bn)}&{\bf F1 (mz)}\\\hline
NER-DE&{\bf 12.50\%}&{\bf 0.685}&0.67&0.583&{\bf 0.727}&0.477&0.012&{\bf 0.601}&{\bf 0.739}&{\bf 0.771}\\
NER-DE (CoreNlp mentions)&13.12\%&0.594&0.584&0.481&0.58&0.466&0.028&0.516&0.663&0.619\\
NER-DE (CogComp mentions)&13.12\%&0.616&0.594&0.484&0.598&{\bf 0.479}&0.006&0.516&0.702&0.656\\
NER-DE (AllenNlp mentions)&35.31\%&0.541&{\bf 0.676}&{\bf 0.606}&0.588&0.335&{\bf 0.032}&0.396&0.675&0.576\\
NER-DE (gold mentions)&11.56\%&0.789&0.817&0.756&0.867&0.488&0.056&0.749&0.815&0.831\\\hline
\end{tabularx}
\caption{Performance of NER-DE system with different methods of mention detection. The left panel shows named people coreference metrics (percentage chains not found and entity F1), and average F1 coreference evaluation combining MUC, \(B^3\) and CEAFE, on all test data. The middle panel shows entity F1 by type of mention, names, pronouns or nominals. The right panel shows entity F1 broken dawn by subgenres in the test data: newswire (nw), broadcast news (bn) and magazines (mz). Existing F1 metric is not sensitive to chains not found, and this explains the difference from the new intuitive NPC metric. For each metric, the best system has been bold-faced}\label{NERDEPeopleMetrics}
\end{table*}

\section{NER-driven Coreference}

The most striking error of current systems is their inability to find a substantial fraction (20\%--40\%) of the NPC entities. In these cases the system may have produced an entity that overlaps with some mentions but does not include the name. Detection of named mentions can be done with high accuracy by named entity recognition systems \cite{DBLP:conf/acl/StoyanovGCR09} and the matching of names can also be done accurately via string matching \cite{Wacholder:1997:DPN:974557.974587,wick2009entity}. Pronoun resolution is more difficult and sometimes even requires background knowledge. However, generally, writing styles dictate how pronouns are used in text. Nominal references, on the other hand, are unlikely to make any noticeable difference since there are only a handful. Inspired by these observations, we present a NER-driven rule based coreference system, NER-DE.  We perform clustering on named entities and link pronouns using basic heuristics. This is simpler than prior work on coreference resolution using clustering. Authors of \cite{cardie1999noun} represent all noun phrases using a set of features and perform clustering on these representations, which can also lead to clusters not having an associated name.

\paragraph{People-mention coreference}
In a first pass, NER-DE finds all spans of text that are {\sc person} named entities, using the Cogcomp NER \cite{RatinovRo09}. We later use a dependency parser \cite{honnibal-johnson:2015:EMNLP} to find the noun phrase of which that named entity serves as the syntactic head. This is the smallest span that contains all descendants of the last word in the named entity. We also use two additional rules: (1) if there is an 'and' immediately after the name, the end of the name is the end of the mention, and (2) if the token immediately before the name is a noun or a proper noun, we include that token and its descendants (in the parse tree) in the mention span. In addition, we experimented with taking mention detection from existing coreference systems rather than using this dependency parsing-based method.

Coreference decisions are done on the named entity spans. We start by focusing on names, initializing a separate chain for each named entity. We then do agglomerative clustering of the names, using a named entity similarity metric~\cite{cogcompnlp2018lrec}. We merge two chains if the longest mention of the first chain has a similarity of more than 0.5 with any mention of the second chain. Most similarity scores are close to 0 and 1 so the system is not sensitive to the chosen similarity cut off.

\paragraph{Resolving pronouns}
We add first-person and third-person singular pronouns to the coreference chains containing named mentions using three rules, restricted by gender compatibility. We consider the gender of the chain to be the gender of the first word in the longest named mention of the chain. We determine the gender using a gazetteer of names, created in an unsupervised manner from Wikipedia. For each person's page, the name is considered female if the number of female pronouns is greater than the number of male pronouns in the first paragraph of the page, else it is considered male. A name is considered unisex if it does not appear in the gazetteer. This gazetteer has been evaluated from names that originate from different languages and gives accuracies in the high 90s. We use the following rules for pronoun resolution:

{\em (i)} If the pronoun is the subject of a verb and the preceding subject is a name, we assign the pronoun to that name. For instance, \textit{he} is assigned to \textit{Henry} in the text - \textit{Henry went to see Barry in the hospital. Afterward, he ate a pizza.},

{\em (ii)} Else, we assign the pronoun to the nearest preceding name such that {\em (a)} the pronoun is not a subject of a verb whose object is the name (or any other name in the same entity cluster) {\em (b)} the pronoun is not an object of a verb whose subject is the name (or any other name in the same chain), and {\em (c)} the name appears at least once in the 100 words preceding the pronoun.

{\em (iii)} If no name satisfies any of the above two conditions, the pronoun is not assigned to a name.

We remove chains that have only one mention, to follow the convention of having only non-singleton chains.

\section{Evaluation of NER-DE}
We evaluate NER-DE on both standard and NPC metrics. The results can be seen in Table ~\ref{NERDEPeopleMetrics}.

NER-DE performs better than the CoreNLP deterministic and statistical systems on the standard metrics. AllenNLP still outperforms all others on this metric. However, as discussed in the earlier sections, this metric is not sensitive to entities not found at all. NER-DE misses 12.5\% of the chains, which is much smaller than what the other systems miss and is still comparatively weaker on the standard metrics. 4.67\% of the chains are missed because NER is not able to find a single named mention in that chain. The rest of the chains are missed due to the coreference algorithm, with no other reference linked to the name.

The existing metrics also do not take into account the resolution of mentions to the correct name, which is incorporated in the NPC F1. NER-DE outperforms all systems on NPC F1, both overall and in the subgenres of nw, bn and mz. It also does better on finding named mentions as well as pronouns. We should note that this improved performance is obtained by basic heuristics of matching pronouns to the nearest preceding named mentions with just a few additional constraints.

We also see that the performance varies by a significant amount depending on the source of mentions used by NER-DE. Surprisingly, the performance of NER-DE is better by using mentions generated using parsing, instead of using mentions generated by CogComp and CoreNLP systems, on both the NPC and standard metrics. The performance of NER-DE using AllenNLP mentions is comparable to parsed mentions using existing metrics but the AllenNLP systems combines mention detection with coreference and thus has better mentions. The NER-DE using parsed mentions still outperforms even the one using AllenNLP mentions on NPC metrics. Using gold mentions gives a very high performance on both the NPC and standard metrics, thus showing that NER-DE, which uses just NER and basic heuristics for pronominal resolution is surprisingly good at named person coreference. Improving mention detection will lead to even better performance.

\begin{table}[b]
\centering
\small
\setlength{\tabcolsep}{4pt}
\begin{tabularx}{0.85\linewidth}{L{0.6}C{0.38}C{0.35}C{0.35}}
\hline
&{\bf Chains not found}&{\bf F1 (NPC)}&{\bf Avg F1 (coref)}\\\hline
PER Allen&31.87\%&0.563&0.71\\
PER NER-DE&12.5\%&0.685&0.67\\\hline
ORG Allen&28.9\%&0.55&0.6\\
ORG NER-DE&28\%&0.41&0.42\\\hline
GPE Allen&13.9\%&0.75&0.75\\
GPE NER-DE&13.2\%&0.63&0.51\\\hline
\end{tabularx}
\caption{Performance in coreference for PER, ORG, and GPE entities}\label{EvalOtherEntities}
\end{table}

\section{Coreference for Other NE types}
Natural language processing applications usually deal with PER, ORG and GPE entities together. Therefore, similar to PER, we evaluated the ORG and GPE coreference chains using standard and NPC metrics. We also built a similar NER-driven rule-based system with agglomerative clustering for resolving names and basic heuristics for resolving pronouns. \ignore{For the ORG and GPE entities, we did not use the dependency parsing-based rules for resolving pronouns but did use dependency parsing for selecting mention spans for named entities.} Unlike PER, the state-of-the-art coreference systems are able to find the same percentage of chains as the NER-driven system and do not suffer from the same issues as person entities that often. They are also able to resolve all types of mentions better than the NER-driven system. These results can be seen in Table \ref{EvalOtherEntities}. \ignore{The percentage of chains not found for the baseline system is a result of (1) NER recall and (2) singleton names, which are removed by the system. In order to determine the effect of each of this, we also computed the percentage of chains in which no entity was found by NER. These percentages are shown in Table \ref{EvalOtherEntities} for each entity. Singletons have a larger effect for each entity type, but the gap in percentage of chains not found between Allen and NER is largest for PER, followed by ORG and GPE.}

\section{Conclusion}
We presented the task of Named Person Coreference (NPC) and showed that the standard coreference metrics are not suitable for the evaluation of this task. We introduced interpretable evaluation metrics that tackle the shortcomings of the standard metrics and also track the different errors made by systems. We showed that the top off-the-shelf systems do not perform well on these metrics. They output many clusters without a link to any name or a link to the incorrect name. We showed that similar issues are not prevalent for other entities. We presented a simple NER-driven algorithm (NER-DE) for the named person coreference task that performs better than top off-the-shelf systems on the new metrics and on par on the existing metrics.


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

As stars coalesce out of dense concentrations within molecular clouds,
ultraviolet radiation, primarily from the most massive
young stars, slowly erodes the clouds away.
Between the ionized region that immediately surrounds young stars and the much colder
ambient molecular cloud, far UV (FUV) stellar radiation with energies between $\sim$ 6~eV and 13.6~eV
dominates the chemistry and heating of the gas.  These FUV photons dissociate molecular hydrogen
in the cloud but do not have enough energy to ionize atomic hydrogen. The resulting
``photodissociation regions" (PDRs) where the cloud absorbs stellar ultraviolet radiation
expand in photoevaporative flows and produce a
wide-variety of pillar shapes and other complex geometries \citep[e.g.][]{hartigan2015}. \par

Newborn stars also inject momentum and energy into molecular clouds
through protostellar outflows \citep[e.g.][]{quillen2005}, stellar winds \citep[e.g.][]{harper2009},
and radiation pressure \citep[e.g.][]{lopez2011,pellegrini2007, lopez2014, krumholz2009,
krumholz2014}. Photoionization and winds are both negative feedback mechanisms in the sense that they
energize and drive turbulence in the parent molecular cloud \citep{matzner2002, krumholz2006, krumholz2014}
and generally reduce the ability of the cloud to create new stars. However, in some cases the
HII region may drive a shock wave into the cloud and compress the gas enough to
trigger new stars to form \citep{bertoldi89,sicilia2014}.
Overall, the feedback mechanisms of photoionization and stellar outflows
explain why star formation rates are much lower than expected from gravitational
collapse alone \citep[e.g.][]{zuckerman1974, krumholz2007, evans2009}. Feedback also
makes up a key component of modern N-body simulations of galaxy
formation \citep[e.g.][]{schaye2015,genel2014,vogelsberger2014,springel2003}. {For a review of feedback processes and their implementation in numerical simulations of star and galaxy formation, see \citet{dale2015}.} \par

In this work, we study the HII regions and PDRs present in three large regions of
massive star formation: Cygnus OB2, NGC 3372 (Carina Nebula) and IC 1396. These
three regions represent massive star formation on very different scales,
allowing us to investigate how different radiation fields from young stars influence molecular clouds.
Cygnus OB2, part of the larger Cygnus X region, is one of the most active regions of massive star
formation in the Milky Way as measured by stellar content. Studies have catalogued at least 169
primary OB stars in Cygnus OB2 with a total stellar mass of $\sim$ $1.6\times 10^{4}$$M_{\odot}$
and {total molecular cloud mass of $\sim$ $7\times 10^5$M$_\odot$ \citep{wright2015,schneider2006,comeron2012}. Due to the high level of extinction, many more early-type stars are likely still undiscovered in the association \citep{berlanas2018}.} Based on its
size, some have even suggested considering Cygnus OB2 as a young galactic globular cluster instead of a
massive star association \citep{knod2000,comeron2002}. Star formation in the Cyg OB2 region appears to have started
$\sim$7 Myr ago and peaked 5 Myr ago \citep{wright2015}, though it
continues in the region today both within isolated globules and as evidenced by a scattered
group of over 1000 protostars \citep{kryukova2014}. { \citet{schneider2016} characterize the population of globules, pillars, and condensations in Cygnus OB2 and suggest an evolutionary path whereby pillars evolve into globules which can evolve into proplyd-like objects. We assume a distance to the Cygnus OB2 region of 1400pc \citep{rygl2012}. }\par

Carina is of a similar scale to Cygnus~OB2, with about
70 O stars catalogued \citep{smith2006}, including the famous LBV star $\eta$ Car. {The cloud mass of the Carina region is roughly $6\times10^{5}$M$_\odot$ with $>10^4$M$_\odot$ at the density sufficient to form more stars \citep{preibisch2012}}. Carina appears to be a
slightly younger 1-3 Myr old region \citep{hur2012} with copious ongoing star formation { that might be triggered by the massive star clusters \citep{smith2010, rocc2013, gac2013}.}
As such, Carina offers an example of a younger analog of Cygnus OB2. {In this paper, we focus on the surroundings of Trumpler 16 and Trumpler 14, the two most massive young clusters \citep{smith2006}. We assume a distance to Carina of 2300pc \citep{smith2006}. }\par

In contrast to Cygnus
OB2 and Carina, IC 1396 has only a few massive stars. Similar to the Orion Nebula,
the primary source of ionization in IC~1396 is a single
O6.5V star (in this case HD 206267, part of the Tr~37 cluster). IC 1396 notably includes
a famous bright rimmed cloud known as the Elephant Trunk Nebula (IC 1396A),
and there is evidence for sequential
star formation in the region. \citet{getman2012} found a spatial age gradient whereby the youngest stars
are located near the cloud, including inside of it, and stellar age increases as one
moves east towards Tr~37. They argue that HD~206267 caused a radiative driven implosion
over the past several Myr that sequentially formed these stars. \citet{sicilia2014} found a very
young protostellar object directly behind the ionization front (IF) consistent with this picture.
Tr~37 itself is found to have an episodic star formation history and is estimated to have a mean age of 4 Myr \citep{sicilia2005, sicilia2015}. {We assume a distance to IC 1396 of 870pc \citep{contreras2002}.}

Two of the brightest emission lines in the near-infrared spectrum of PDR's,
H$_{2}$ 1-0 S(1) and Br$\gamma$, have proved to be excellent tracers of PDR interface shapes
because they highlight where the FUV radiation is absorbed while penetrating through most of
the obscuring dust along the line of sight.  The 2.12$\mu$m H$_{2}$ emission line
is a fluorescent line of molecular hydrogen that is emitted in a cascade following the
absorption of a FUV photon.  Photons that have energy above roughly 10eV can excite an
electron in H$_{2}$ into the Lyman and Werner bands, and, as the electron cascades
to the ground state, there is a chance it passes through the transition that emits
a 2.12$\mu$m photon.  Hence, this line traces exactly where FUV radiation is absorbed in the cloud.
In contrast, the 2.17$\mu$m Br$\gamma$ line is a
recombination line of ionized hydrogen and thus traces the H$^{+}$ in the regions.
Because the lines are bright and differential extinction is negligible,
once any K-band continuum is removed the resulting images in Br$\gamma$ and H$_2$
allow us to probe the physical structure of the ionization fronts and PDRs throughout
a region of star formation.  These emission lines will undoubtedly be primary choices
for this type of work once JWST is launched.

In this paper we present extensive maps of
continuum-corrected H$_{2}$ 1-0 S(1) and Br$\gamma$ line emission in Cygnus OB2, Carina, and
IC 1396. Motivated by the distinct spatial offset between the H$_{2}$ and Br$\gamma$
emission observed in Carina's PDR interfaces \citep{hartigan2015}, we identify roughly 450
locations along the PDR interfaces in the three star forming regions that are suitable for
extracting high signal-to-noise spatial profiles in both lines.
We then use numerical PDR models to understand how the size of the offset is determined by
the geometries and radiation fields present in these regions.
In Section \ref{sec:data} we discuss the data acquisition
and reduction; in Section \ref{sec:Analysis}, we discuss the regions and the offsets we find in
each; and in Section \ref{sec:models} we introduce the numerical models and parameters used to compare with
and interpret our observations. Section \ref{sec:conclusions}
 summarizes our results and conclusions.

\begin{deluxetable}{lcc}
\tablecaption{Summary of Observations\label{tab:observations}
}
\tablehead{
\colhead{Region} & \colhead{Dates} & \colhead{Telescope}
}
\startdata
Cyg OB2 & Sept 19-22, 2008 & 4m Mayall, KPNO\\
           & Nov 12-20, 2012 & 4m Mayall, KPNO\\
Carina     & March 11-18, 2011     & 4m Blanco, CTIO\\
IC 1396    & Oct 27-30, 2012     & 4m Mayall, KPNO\\
           & Nov 12-20, 2012     & 4m Mayall, KPNO\\
\enddata
\end{deluxetable}

\section{Data} \label{sec:data}

\begin{figure*}
\includegraphics[width=0.95\textwidth]{car.pdf}
\caption{H$_{2}$ mosaic of Carina region from \citet{hartigan2015}. Irradiated molecular interfaces appear as the dark pillars
and walls. The continuum has not been subtracted. The WR, O, and early
B stars from \citet{gagne2011} are shown in the black circles. {The molecular cloud interfaces analyzed in this work are indicated by the red boxes}. North is up and East is to the left in all the images. }
\label{fig:carina_img}
\end{figure*}

\subsection{Reduction}

The near infrared (NIR) data used in this project were all taken with the NOAO
NEWFIRM instrument. The imager consists of four 2048x2048 pixel detectors with a gap of 35
arcseconds between chips \citep{probst2008}. Attached to the NOAO 4m telescopes, the imager has a pixel
scale of 0.4 arcsec/pix. Details for the data acquisition for the three regions are given
in Table \ref{tab:observations}. NEWFIRM has a data
reduction pipeline \citep{swaters2009} that attempts a sky subtraction by using a median filter over the
nearest 3 images in time before and after a frame is taken. As discussed in
\citet{hartigan2015}, the pipeline is inadequate for regions of large-scale nebulosity.  We therefore re-reduce the data
in the manner outlined in \citet{hartigan2015}. To improve the S/N of our K-band mosaic, we used
an archival NEWFIRM image of IC 1396A taken in 2009 by T. Megeath.
This improved the S/N of the tip of IC 1396A significantly but did not help
for regions further down the trunk. The maps of Carina were presented and thoroughly discussed
in \citet{hartigan2015}, and several of the H$_2$ images of globules in Cygnus OB2 were shown
in \citet{hartigan2012}.

\begin{figure*}
\includegraphics[width=0.95\textwidth]{cyg.pdf}
\caption{H$_{2}$ mosaic of Cygnus OB2 region. The location of the WR, O, and early B
stars from \citet{wright2015} are shown in the black circles. {The molecular cloud interfaces analyzed in this work are indicated by the red boxes}. The continuum has not been subtracted.}
\label{fig:cyg_img}
\end{figure*}

The final stacked results were then aligned between the three filters. S/N and exposure time
varied within the frames due to the complicated dither patterns. The Carina data had average
seeing of 1.09", IC 1396 had 1.1", and Cyg OB2 had 1.1".
The images were acquired through a variety of sky conditions and are not flux-calibrated, but
this aspect does not affect the spatial offsets observed between H$_{2}$ and Br$\gamma$ emission
or the spatial profile of either line, which are the focus of this paper.

\subsection{Continuum Subtraction}

The H$_{2}$ and Br$\gamma$ images are contaminated with continuum emission
that arises primarily from starlight reflected from dust in the molecular clouds.
We can remove the continuum by subtracting an appropriately scaled K-band image from
both the Br$\gamma$ and H$_2$ mosaics. Because the wavelengths of H$_2$ (2.12$\mu$m) and Br$\gamma$ (2.17$\mu$m)
are nearly identical and the continuum arises from the Rayleigh-Jeans part of the black-body
spectrum for a usual star, the amount of continuum is nearly proportional to the bandpasses of
each narrow-band filter. In this situation, there are three unknowns: the isolated
H$_{2}$ and Br$\gamma$ emission and the pure continuum emission, and three constraints: the total emission
recorded in each of the three filters. It is therefore possible to do filter
``algebra'' to isolate the line emission. In particular, we can write:

\begin{equation}
H_{2}^{tot} = H_{2} + \alpha K \\
\label{eq:h2}
\end{equation}

\begin{equation}
Br\gamma^{tot} = Br\gamma + \beta K \\
\label{eq:brg}
\end{equation}

\begin{equation}
K^{tot} = H_{2} + Br\gamma+  K \\
\label{eq:k}
\end{equation}

\noindent
where X$^{tot}$ is the observed flux in the X-filter image and X is the contribution of the
emission line. The constants $\alpha$ and $\beta$ are related through the bandpass widths of the H$_{2}$
and Br$\gamma$ filters, respectively. We determine their ratio empirically by constructing
scaled subtractions between
the H$_{2}^{tot}$ and Br$\gamma^{tot}$ images. Stars will generally not have any line emission in H$_{2}$
or Br$\gamma$ and will subtract out with the correct values of $\alpha$ and $\beta$.
We find that H$_{2}^{tot}$ - $C\times$Br$\gamma^{tot}$, where $C\equiv\alpha/\beta$, optimally removes stars when $C$=1.19 for all regions in
Carina and Cygnus except for one grid in Cygnus where $C$=1.40 is required. For IC 1396,
$C$ was found to be between 0.85 and 1.05 depending on the region of the
mosaic. The scaling factor includes the effects of different
filter bandpasses, exposure times, and sensor sensitivities.

For some constant, $D$, adding $D$ times Eq. \ref{eq:brg} to Eq. \ref{eq:h2} and subtracting $D$
times Eq. \ref{eq:k} yields:

\begin{equation}
H_{2}^{tot} + DBr\gamma^{tot} - DK^{tot}= H_{2}(1-D) + K(\alpha + \frac{\alpha}{C}D-D)
\label{eq:deq}
\end{equation}

When $D$ is chosen such that the stars are properly subtracted, we see that the
coefficient of $K$ on the right hand side of Eq. \ref{eq:deq} must be zero, implying:

\begin{equation}
\alpha = \frac{D}{\frac{D}{C}+1}
\end{equation}

$\beta$ can then be found from the measured values of $C$. Once we know $\alpha$ and $\beta$, the line emission is given by:

\begin{equation}
H_{2} = \frac{H_{2}^{tot}-\frac{\alpha}{1-\beta}\times(K^{tot}-Br\gamma^{tot})}{1-\frac{\alpha}{1-\beta}}
\end{equation}

\begin{equation}
K = \frac{K^{tot} - H_{2} - Br\gamma^{tot}}{1-\beta}
\end{equation}

\begin{equation}
Br\gamma = Br\gamma^{tot} - \beta K \\
\end{equation}

$\alpha$ and $\beta$ are slightly different for the different mosaics but always close to the
ratios in bandpass between the H$_{2}$/Br$\gamma$ filters and the K filter. Our $\alpha$ and
$\beta$ for Carina correspond very closely to those found by \citet{yeh2015} who imaged 30 Dor
with NEWFIRM.

\begin{figure}
\includegraphics[width=0.45\textwidth]{ic1396_img.pdf}
\caption{H$_{2}$ mosaic of IC 1396 region. Location of HD206267 in the center of Trumpler 37
is denoted by a black circle. IC 1396A, the Elephant Trunk nebula, is visible across the bottom of the frame pointing towards HD206267. Part of the edge of the cavity is visible in the upper right. {The molecular cloud interfaces analyzed in this work are indicated by the red boxes.} The continuum has not been subtracted.}
\label{fig:ic_img}
\end{figure}

\section{Analysis} \label{sec:Analysis}
{Our primary goal in this paper is to quantify how ionization fronts (IFs) and
photodissociation regions (PDRs) appear in Carina, Cyg OB2, and IC 1396.
In this section we summarize the general characteristics of the H$_2$ emission-line
nebulae in each region (Section \ref{sec:h2_overview}), investigate how the spatial
profiles of Br$\gamma$ emission behave near the cloud surfaces (Section \ref{sec:brgProfiles}) and
consider the observed spatial offsets between H$_2$ and Br$\gamma$ across the molecular
cloud surfaces (Section \ref{sec:offsets}). Our study quantifies differences in the size of this
offset between the three star forming regions as well as trends in offset size
with cloud geometry and irradiation environment.}

\subsection{The PDR Interfaces in Carina, Cyg OB2, and IC 1396}\label{sec:h2_overview}

{The H$_2$ images are excellent tracers of PDR interfaces, and reveal
a variety of pillars, walls, and globules in Carina, Cyg OB2, and
IC 1396 (Figs. \ref{fig:carina_img}-\ref{fig:ic_img}). Each star-forming region has its own general
character in the H$_2$ images. The image of Carina
in Figure \ref{fig:carina_img} contains several clusters, including Tr 14, the bright
LBV star $\eta$ Car, and a wide variety of pillars, walls, and globules.
In contrast, though Cyg~OB2 also has large-scale H$_2$ walls, the brightest
PDRs are associated with irradiated globules located in the vicinity of the
concentration of massive stars in the middle of the images.}

{The last region, IC 1396A, is dominated by the Elephant Trunk
Nebula, a long pillar structure.} Also present are the edges of a large, circular bubble surrounding Tr 37.
Because of the low S/N of
the K-band image, we focus primarily on the PDR interfaces along and at the base
of the trunk where the S/N is the best in the K-band. In all of our regions
the Br$\gamma$ emission is located closer to the ionizing sources than the H$_{2}$,
compatible with the picture of molecular clouds being ionized and dissociated by radiation from the
newborn stars. For many of the bright interfaces, streams of photoionized gas are visible
in Br$\gamma$ as they evaporate from the surface of the clouds.

\begin{figure*}
\includegraphics[width=0.99\textwidth]{median_profile_shapes.pdf}
\caption{{Continuum subtracted Br$\gamma$ profiles for five PDR interfaces in Carina.
The left column shows all Br$\gamma$ profiles for a given interface. The middle column shows the
median profile shape along with the best power-law fit. The right column displays an
H$_2$ image of the object (not continuum-subtracted) together with the locations of all the
slices (shown in blue) that were used to fit the Br$\gamma$ profiles. The labels of the interfaces, as listed in the appendix, are indicated in the upper left.}}
\label{fig:prof_shapes}
\end{figure*}

\begin{figure}
\includegraphics[width=0.46\textwidth]{powerlaw_slopes.pdf}
\caption{Histogram of power law indices for 26 different interfaces in Carina and Cyg OB2, broken
into groups by their approximate geometric shape. }
\label{fig:powerlaw}
\end{figure}

\subsection{Br$\gamma$ Profiles and Interface Shapes}
\label{sec:brgProfiles}

We can study the IFs and photoevaporative flows in these regions
by measuring the spatial profiles of Br$\gamma$ intensity near the surface of the clouds.
Br$\gamma$ emission arises from photoevaporative flows off of the cloud.
Its brightness should reflect the density profile of the flow because Br$\gamma$ emissivity $\propto
n_e n_p$. Photoevaporative flows ablate from cloud surfaces at a few
times the sound speed \citep{bertoldi89}. If we assume that the flow moves with a roughly
constant velocity then conservation of mass requires that the density falls like $n\propto r^{-2}$ for
spherical clouds and $n\propto r^{-1}$ for cylindrical clouds, where $r$ is the distance from the
surface of the cloud. Assuming that the transition from neutral to ionized material is over a short distance, we then expect the Br$\gamma$ emissivity to follow a power law profile of the form:
\begin{equation}
\epsilon_{\rm Br\gamma} \propto \left(\frac{r+r_c}{r_c}\right)^{-b}
\label{eq:brg_emissivity}
\end{equation}

\noindent
where $r$ is the distance along the profile with $r=0$ corresponding to the peak in
Br$\gamma$ emission at the cloud's surface, $r_c$ is the radius of curvature of the cloud, and $b=2$ for a pillar and $4$ for a globule. Assuming that the photoevaporative flow extends outward from the cloud's surface to infinity, projecting the Br$\gamma$ emissivity (Equation \ref{eq:brg_emissivity}) onto the plane of the sky yields an observed Br$\gamma$ intensity profile of the form:
\begin{equation}
I_{\mathrm{Br\gamma}} = I_{\mathrm{max}}\left(\frac{r+r_c}{r_c}\right)^{-a}
\label{eq:brgProfile}
\end{equation}

\noindent
where $I_{\rm max}$ is
the peak Br$\gamma$ intensity on the profile and $a=b-1$. In other words, the effect of projection on the Br$\gamma$ profile is to make the power law shallower, decreasing the power law exponent by one, regardless of whether it is a pillar or globule. In the more realistic case that the photoevaporative flow does not extend to infinity and a uniform background takes over at some distance from the cloud, the projected intensity profile will not strictly be a power law since it goes to a constant at some radius. Assuming that the flow extends out at least a few times the radius of curvature of the cloud from the cloud's surface, the profile can be approximated as a power law with the exponent in between the unprojected and infinite flow cases. This means that we expect globules to have roughly power law Br$\gamma$ intensity profiles with exponent in between $-3$ and $-4$ and pillars to have exponents between $-1$ and $-2$\footnote{For instance, numerically projecting Equation \ref{eq:brg_emissivity} onto the plane of the sky for a flow that extends to $2r_c$ from the cloud's surface yields a projected Br$\gamma$ intensity profile of $I_{\rm Br\gamma} \propto r^{-1.78}$ for a pillar and $I_{\rm Br\gamma} \propto r^{-3.22}$ for a globule}.\par

To compare this with our data, we extract slices from the continuum-subtracted Br$\gamma$ images across cloud interfaces.
Slices were taken only on interfaces that appear edge-on, and we averaged the profiles over 5
pixels perpendicular to the slice to improve S/N. Slices must exclude stars,
and this requirement limits the total number of slices.
Slices were taken on a variety of interfaces, including pillars, walls, and globules.
For a given interface, we extracted several independent slices and
aligned them by using the peak Br$\gamma$ emission. These slices generate a median profile
shape for each interface. Only the part of the slices on the HII region side of the Br$\gamma$ peak is kept for this analysis. To isolate the Br$\gamma$ emission in the photoevaporative flow from foreground/background nebulosity, a background level is subtracted from each slice.
The background is quantified as the minimum pixel value along the slice. To make sure that the slices extend far enough out into the HII region that this minimum value does represent a background level and is not due to the photoevaporative flow itself, we further restrict slices to be longer than the radius of curvature of the object.
 We grouped interfaces into three broad categories: globules, which have roughly spherical symmetry; pillars, which
have approximately cylindrical symmetry; and intermediate shaped clouds, which are a mixture of the
two.  The pillars are specifically defined as those objects with a long axis being at least twice the length of a short axis.
 {For all the interfaces we can estimate a radius of curvature directly from the H$_2$ images by fitting an arc with
a constant radius of curvature to the feature.} The profiles were fit in the range $r=0$ to $r=r_c$.
To account for the blurring effect of seeing, the profile of Equation \ref{eq:brgProfile} is first convolved with a Gaussian of the width of the seeing before being fit to the data.\par

Figure \ref{fig:prof_shapes} shows these fits for five typical interfaces in Carina. The middle column shows the median
profile shape and the best fitting power-law. Images of the interfaces are also shown. The
top two interfaces are categorized as `pillars', the middle one as `intermediate', and the bottom
two as `globules'. The power law slope is steeper for the globules than the pillars,
as expected, because the extra curvature in the globules causes the gas to diverge over
shorter spatial scales than in the pillars. Doing this analysis for 26 objects in Carina and
Cygnus OB2\footnote{The interfaces in IC 1396 had large radii of curvature which made it impossible to get slices that were as long as the radius of curvature while still avoiding stars.} we derive the power law fits shown in Figure \ref{fig:powerlaw}. The
`globule' interfaces have an average power law index of $-3.45\pm0.22$ whereas the `pillar' interfaces have
an average index of $-2.54\pm0.16$. A two sample $t$-Test indicates that these means are significantly different with a $p$-value of 0.004. The globule power law exponents agree well with the expectation above that they fall in between $-3$ and $-4$, and the globule objects have steeper profiles than the pillars, also
consistent with simple theoretical expectations. The power law exponents for the pillars are likely greater than the expected values between $-1$ and $-2$ because the objects are not exclusively cylindrical but consist of a mixture of cylindrical and spherical curvature.  A list of all the interfaces considered can be found in Appendix \ref{sec:appendix}

\begin{figure*}
\includegraphics[width=0.99\textwidth]{sample_slice_ic.pdf}
\caption{{H$_2$ image (not corrected for continuum) of one of the objects in IC 1396
used for measuring spatial offsets between H$_2$ and Br$\gamma$ emission.
The green lines mark locations where the continuum-corrected
spatial profiles were extracted for the emission lines. The spatial
offset between H$_2$ and  Br$\gamma$ is shown at right for the position
defined by the purple line. The interface's label, as given in the appendix, is indicated in the upper left.}
\label{fig:slice_ic}}
\label{example_slice_ic}
\end{figure*}

\begin{figure*}
\includegraphics[width=0.99\textwidth]{sample_slice_car.pdf}
\caption{{Same as Figure \ref{fig:slice_ic} for an object in Carina.}
\label{fig:slice_car}}
\end{figure*}

\begin{figure*}
\includegraphics[width=0.99\textwidth]{sample_slice_cyg.pdf}
\caption{{Same as Figure \ref{fig:slice_ic} for an object in Cygnus OB2.}
\label{fig:slice_cyg}}
\end{figure*}

\subsection{Spatial Offsets Between H$_2$ and Br$\gamma$}
\label{sec:offsets}
As noted in \citet{hartigan2015}, there is a distinct offset between the H$_{2}$ and Br$\gamma$ emission
when slices are extracted across PDR interfaces. Since the H$_{2}$ and Br$\gamma$ lines are so
close in wavelength, this offset should be independent of extinction, and it should reflect the physical
size of the PDR and ionization front. In this section we
quantify this offset in Carina, IC~1396, and Cyg~OB2. Sample slices for each of the three regions are shown in Figs. \ref{fig:slice_ic}-\ref{fig:slice_cyg}, displaying typical offsets between H$_{2}$ and Br$\gamma$ for irradiated clouds.

Because the three star formation regions we are analyzing
are at different distances, the images will have markedly
different spatial resolutions even if the seeing were the same.
To correct for this effect, we bin the IC 1396 and Cygnus OB2 data so that they have the same
pixel scale as the Carina data (Carina is the furthest region of the three).
We assume a distance of 1400~pc to Cyg OB2 \citep{rygl2012},  2300~pc to Carina \citep{smith2006}, and
870~pc to IC 1396 \citep{contreras2002}. At the distance of Carina, one pixel corresponds to $1.38\times10^{16}$cm. To equalize the
seeing, we blurred all the data down to the level of the image with the
worst seeing by convolving the images with Gaussians of varying widths. The image with the
worst seeing is the K-band Carina image with a FWHM seeing of $3.95\times10^{16}$cm (at the
distance of Carina), and all final images have this same level of smoothing.

To quantify the observations we simply measure the spatial offset between
the peak emission in Br$\gamma$ and in H$_2$ in the smoothed images.
Since the slices do not need to be longer than the interface's radius of curvature, we can consider more
distinct objects than in Section \ref{sec:brgProfiles}. In total, we consider 33 objects across the three regions consisting of around 450 distinct slices. These objects are listed in Appendix \ref{sec:appendix}.\par

Figure \ref{fig:offset_3regions} shows the
cumulative distribution for physical offset sizes for the three regions.
IC 1396 appears to have the largest offsets while Cyg OB2 and Carina
have roughly similar offsets. IC 1396 has median offset size of $9.9\times10^{16}$cm, while Carina has median offset size of $7.9\times10^{16}$cm, and Cyg OB2 has median offset of $7.8\times10^{16}$cm.
Figure \ref{fig:offset_3regions} shows that the offsets are all within roughly a factor of two of
$10^{17}$cm, despite the fact that the three regions span at least
three orders of magnitude in the intensity of ionizing radiation.

{We can use our data to determine if the observed offsets depend
upon the radius of curvature of the cloud or the strength of the radiation field.
Many interfaces have approximately spherical or pillar-shaped morphologies,
while others appear intermediate between these shapes. The H$_2$ image gives a characteristic radius of curvature for each interface, with a typical error
of $\pm$ 10\%. }To estimate the intensity of ionizing radiation incident onto each interface,
we use catalogs of the O and early B-type stars in the regions. For Carina,
we reference the catalog of O/B stars compiled by \citet{gagne2011}, and for Cygnus OB2, we apply
the list of \citet{wright2015}. Because HD~206267 is the lone O-star in IC 1396, we need
only consider that source.  We convert the spectral type to ionizing photon output
using values from \citet{martins2005}. To approximate the flux of ionizing photons incident at
each slice location, we consider only the O/B stars that are visible to that slice
in the direction above the cloud interface.
We combine the stars via:

\begin{equation}
F = \sum_i \frac{Q_i \cos(\theta)}{4\pi r^2},
\end{equation}

\begin{figure}
\includegraphics[width=0.45\textwidth]{offsets_global.pdf}
\caption{Distribution of offset sizes.
\label{fig:offset_3regions}}
\end{figure}

\begin{figure}
\includegraphics[width=0.45\textwidth]{off_v_f.pdf}
\caption{Measured offset sizes versus the estimated incident flux of H-ionizing photons.
Objects in Carina are shown in red, those from Cygnus OB2 are in blue, and those from IC 1396
are in green. The flux of ionizing photons spans three orders of magnitude.}
\label{fig:o_v_f}
\end{figure}

\begin{figure}
\includegraphics[width=0.45\textwidth]{off_v_rc.pdf}
\caption{Measured offset sizes versus the measured radius of curvature. Objects in Carina are shown in
red, those from Cygnus OB2 are in blue, and those from IC 1396 are in green.}
\label{fig:o_v_rc}
\end{figure}

\begin{figure}
\includegraphics[width=0.46\textwidth]{plane_fit.pdf}
\caption{Edge-on view of best fitting plane to the measured offset values along with best-fitting parameter
values. Objects in Carina are shown in red, those from Cygnus OB2 are in blue,
and those from IC 1396 are in green.}
\label{fig:plane_fit}
\end{figure}

\noindent
where $Q_i$ is the ionizing photon output of the $i$-th star, $\theta$ is the angle
between the line connecting the star to the interface and the normal to the interface
(which is taken to be the slice we make), and $r$ is the projected distance
between the star and the interface.

To account for projection effects, we calculate an uncertainty in
the flux, $F$, by assuming that the distance, $r$, is known
to within 25\%\footnote{Which corresponds to the vector between the star and interface being within 37$^{\circ}$ to the plane of the sky.}. This uncertainty is then propagated through to the final flux value
for each slice. We average the flux values for each slice on an object to get an estimate for the flux value for that object. Since the uncertainties for the flux value at each slice on an object are not uncorrelated, we just take the rms of these uncertainties as an estimate for the uncertainty in the ionizing flux for that object. { It is also possible to estimate the FUV flux at each point in the regions from the observed far-infrared (FIR) intensities \citep[e.g.][]{kramer2008, rocc2013} but since the EUV flux is more relevant to studying the IFs and photoevaporative flows, we do not explore that method here.}

The offsets measured from slices
on an object are averaged to give a mean offset for that object and a
standard error of the mean.

After these measurements are completed, each object has an associated radius of
curvature, an ionizing flux level, and a spatial offset between Br$\gamma$ and H$_2$.
Figure \ref{fig:o_v_f} shows the measured offsets as a function of the estimated
incident flux. The range in $\log(F)$ shows that our objects differ in the level of
incident flux by around three orders of magnitude between the weakest irradiation environment in IC
1396 and the strongest in Carina. Figure \ref{fig:o_v_rc} shows the offsets as a function of
the measured radius of curvature. There is a reasonably strong positive correlation between
the offset size and the radius of curvature and a very weak negative correlation with respect to
the incident flux. A Spearman rank-order test for Figure \ref{fig:o_v_f} gives a $p$-value of no correlation of 0.05 and a $p$-value of 0.007 for Figure \ref{fig:o_v_rc}. The low $p$-value in Figure \ref{fig:o_v_f} is probably driven by the IC 1396 interfaces at low incident flux. Not including the IC 1396 interfaces yields $p$-values of 0.3 and 0.07 for Figures \ref{fig:o_v_f} and \ref{fig:o_v_rc}, respectively.
 Both Figures \ref{fig:o_v_f} and \ref{fig:o_v_rc} show a significant spread in the
offset sizes even after all the individual slices for an object have been
averaged. Each point in Figures \ref{fig:o_v_f} and \ref{fig:o_v_rc} is already an
average of 10-20 individual slices.

The spatial offsets between Br$\gamma$ and H$_2$ should depend jointly upon both the incident
flux and radius of curvature, which introduces scatter to Figures \ref{fig:o_v_f} and
\ref{fig:o_v_rc} when offsets are plotted against either $\log(F)$ or
$\log(r_c$) individually. To correct for this effect, we fit a plane to the offset values of the form:

\begin{equation}
\mathrm{O} = a\log(F) + b\log(r_c) + c
\end{equation}

\noindent
where $\mathrm{O}$ is the measured offset. Figure \ref{fig:plane_fit} shows an edge-on view of this plane
along with the best fit parameters $a,\; b,$ and $c$.  The scatter in the planar view
(Figure~\ref{fig:plane_fit}) is less than that for $\log(F)$ (Figure~\ref{fig:o_v_f}) and
for $\log(r_c)$ (Figure~\ref{fig:o_v_rc}), though considerable intrinsic scatter remains,
caused in large part by the complex geometrical shapes of the individual objects and uncertainties in the true (deprojected)
distances of the ionizing sources from the cloud interfaces.
Nevertheless, there is a weak positive trend, $b = 2.08\pm 0.44$, of the offsets with radius of curvature
of the cloud and a weak negative trend for the offset with the incident flux,
$a = -1.47\pm0.21$.

\section{Photoionization Models}\label{sec:models}

\subsection{Motivation}

{Our two primary motivations for developing a numerical model of these regions are (1) to explain the
differences and similarities in the offset sizes observed in the three regions and (2) to explore the relevant
physical processes that define the size of the offsets.} The fact that the three regions show
similar offsets indicates that offsets might end up being more useful as a rough distance
indicator than as a PDR diagnostic tool.
Hence, rather than use the spatial offset between Br$\gamma$ and H$_2$
to determine physical conditions, the focus will be to understand
the underlying processes that create these offsets and to predict any correlations with other observables.

\begin{figure*}
\centering
\includegraphics[width=.92\textwidth]{globule2.pdf}
\caption{Schematic of a photoevaporative flow from an irradiated globule.
Top: Ionizing flux $F_{\rm inc}$ from the left impinges upon a dense spherical globule
of radius $r_c$ and drives an outflow. The green colors depict
Br$\gamma$ emission and the red show H$_2$ emission. Bottom: Recombining gas in the
outflow reduces the incident ionization flux by a factor of $q$ at the ionization front.
The ionization front progresses into the cloud at a velocity $v_{\rm IF}$ $\sim$ $v_{D}$,
where $v_{D}$ is the velocity of the isothermal D-shock, set by the
sound speeds $c_i$ and $c_n$ in the ionized and neutral gases, respectively.
The model solves ionization and dynamical equations
for $q$ to determine $n_{\rm I}$, which in turn becomes an input parameter to
\textsc{Cloudy} for determination of the H$_2$ emission. The velocity arrows
depicted in the figure are relative to the unperturbed globule.}
\label{fig:model_setup}
\end{figure*}

\begin{figure*}
\includegraphics[width=0.95\textwidth]{model_x_n.pdf}
\caption{Model results for models with $r_c = 5\times10^{17}$cm. The left plot shows profiles of
the ionization fraction. The right plot shows the electron number
density for the same models. In order of increasing ionization fraction at a given radius,
the models have $\log(F) = 9.3, 9.8, 10.3, 11.3, 12.3, 13.3$ cm$^{-2}$s$^{-1}$. }
\label{fig:model_1}
\end{figure*}

\subsection{Model of a Photoevaporative Flow}

The Br$\gamma$ emission originating near the molecular clouds comes from a photoevaporative flow
of ionized gas off the cloud surface. Therefore, to accurately model the offsets and their
dependence on various parameters, a rough model of a photoevaporative flow is necessary.
We follow \citet{bertoldi89} and solve the equations of hydrodynamics and ionization balance in a simplified
spherically symmetric, 1-D geometry and in the steady-state. In this formalism, all partial time derivatives
in the fluid equations are set to zero. A picture of the physical situation along with a definition of the variables is shown in Figure \ref{fig:model_setup}. We solve the continuity equation, momentum equation,
and equations for ionization balance as:

\begin{equation}\begin{aligned}
\frac{d}{dr}\left(r^2n_{\rm H}v\right) &= 0 \\
v\frac{dv}{dr} &= -\frac{1}{\rho}\frac{dP}{dr} \\
\frac{dF}{dr} &= (1-x) n_{\rm H}\sigma_\nu F \\
\frac{dx}{dr} &= \frac{1}{v}\left(-\alpha_b x^2 n_{\rm H} + (1-x)\sigma_{\nu} F\right)\\
\end{aligned}\end{equation}

\noindent
where $n_{\rm H}$ is the hydrogen number density, $\rho = m_{\rm H} n_{\rm H}$, $P$
is the pressure, $v$ is the velocity of the flow, $F$ is the flux of
ionizing photons at a certain point in the flow, $x$ is the ionization fraction, $\alpha_b$
is the case-B recombination coefficient, and $\sigma_{\nu}$ is the frequency averaged
photoionization cross section for hydrogen. We use $\alpha_b=2.7\times10^{13}$cm$^{3}$s$^{-1}$
and $\sigma_{\nu} = 3.0\times10^{-18}$cm$^{2}$. Following \citet{lefloch1994},
instead of solving a separate energy differential equation, we write the sound speed
as:

\begin{equation}\begin{aligned}
c^{-2} = c_i^{-2} + \left(c_n^{-2} - c_i^{-2}\right)\left(1-\sqrt{x}\right)^2
\end{aligned}\end{equation}

\noindent
which is essentially an arbitrary interpolation function between the two extremes, $c_n$ and $c_i$,
the isothermal sound speeds for the cold neutral gas and the warm ionized gas,
respectively. We take $c_n$ = 0.8 km/s corresponding to $T=100$~K neutral hydrogen and $c_i$=11.4 km/s
corresponding to $T=10^{4}$~K ionized hydrogen. As discussed in \citet{lefloch1994},
writing the sound speed explicitly as such requires that the flow be in thermal equilibrium. More specifically, it requires that
the thermal timescale be much less than the dynamical timescale of the flow. \citet{lefloch1994} argue
that this is the case for flows with strong irradiance where the flow acts as
an insulating boundary layer and absorbs most of the incoming ionizing photons before they
reach the IF. This will be the case for the regions we consider with strong irradiance in Carina and Cygnus but
might not be accurate for the objects in IC 1396 which have fairly weak incident radiation.

We assume the ionization front is D-type and is therefore preceded by a shock wave
in the steady-state (Figure \ref{fig:model_setup}). Following \citet{bertoldi89} we model the flow from the neutral part outwards.
\citet{bertoldi89} found that the velocity of the IF
is approximately the D-critical velocity under a
wide range of conditions. We write the IF velocity as

\begin{equation}\begin{aligned}
v_{\rm IF} = \phi_d v_D
\end{aligned}\end{equation}

\noindent
where $v_D$ is the D-critical velocity $\approx c_{n}^2/2c_i \approx 0.03$ km/s \citep{spitzer1968} and $\phi_d$ is an order unity number.
We will assume that $\phi_d = 1$. Additionally,
the density of the shocked gas layer will be

\begin{equation}\begin{aligned}
n_{\rm I} = \frac{F_{\rm IF}}{v_{\rm IF}}
\end{aligned}\end{equation}

\noindent
where $F_{\rm IF}$ is the flux of ionizing photons that actually reaches the IF. We
parameterize the effect of shielding by the flow by writing $F_{\rm IF} = F_{\rm inc}/q$
where $q$ is unknown \textit{a priori} and $F_{\rm inc}$ is the incident flux. The model
starts in the neutral part of the cloud with $F_0 = 10^{-8}F_{\rm inc}$ and $x_0 = 10^{-8}$
and numerically integrates the equations as it moves outwards. The gas accelerates as it
heats up and becomes super-sonic. However, since we are interested in the density profiles near
the IF and not further out along the flow, we do not follow the super-sonic
evolution. Once the gas becomes faster than $0.99$ times the sound speed, we set the gas velocity as the
sound speed and do not consider the momentum equation any further. This allows us
to avoid the complication that the momentum equation is actually singular at the sonic point.
The models are integrated from $1r_c$ to $5r_c$. The parameter $q$ is chosen so that
the flux at $5r_c$ is within 5\% of the intended incident flux, $F_{\rm inc}$. We
repeat the process for various radii of curvatures and fluxes relevant to our regions. Ionization fraction
and density profiles for models with $r_c = 5\times10^{17}$cm for a variety of incident ionizing fluxes are shown in Figure \ref{fig:model_1}. \par

The models show a spike in the density of HII in the profiles in the ionization fronts.
This occurs because there is a large density contrast across the IF of roughly a
factor of $c_i/v_{\rm IF} \approx 200$. This causes the IF to have higher density than
the photoevaporative flow and so even if the IF is not completely ionized, it can
have a higher HII density than in the flow. {This spike is model dependent (it is sensitive to our simplified energy equation) and it is unclear if it would be present if heating and cooling are explicitly considered. We do not consider this here as the spike is too narrow to affect the observed offset between Br$\gamma$ and H$_2$. It will be completely smoothed once blurred by the seeing.} \par

Figure \ref{fig:model_2} shows the scalings of
the maximum HII density in the flow with the flux and radius of curvature. Also
shown is the expectation for a flow with constant velocity and infinitely thin IF given
by \citep[e.g.][]{tielens05}:

\begin{equation}\begin{aligned}
n_{\mathrm {HII,\; max}} = \left(\frac{3 F}{r_c \alpha_b}\right)^{1/2}
\\
\label{eq:ana}
\end{aligned}\end{equation}

This density is found by assuming that the photoevaporative flow acts as an insulating boundary
layer that absorbs most of the ionizing photons before they reach the IF. We
can see that the scalings in the simulations with regard to $F$ and $r_c$ are
the same as Equation \ref{eq:ana} but that the simulation densities are roughly $1.5\times$ higher. This
is due to the fact that the ionization fronts have non-zero thickness. Equation \ref{eq:ana} assumes
that all of the ionized hydrogen lies outside the IF and moves at
a constant velocity where, {in our models}, the peak density of HII occurs inside the IF where
gas is not completely ionized but is very dense because the gas is still moving
slowly. As mentioned above, conservation of mass requires that the density of the gas before
the IF is roughly $\approx 200$ times that of the gas on the ionized side
of the IF. Within the IF, the gas density will be somewhere in between. \par

This scaling of $n_{\mathrm {HII,\; max}}\propto F^{1/2}$ is found to agree quite well with the pillar mass loss rates in Carina, M16, and NGC 3603 analyzed in \citet{mcleod2016}.

\begin{figure}
\includegraphics[width=0.48\textwidth]{model_scaling.pdf}
\caption{Solid lines show the maximum HII density in the simulations as a function of $F$
and $r_c$. In dotted lines is shown the expectation for a photoevaporative flow with constant
velocity and infinitely thin IF (cf. Equation \ref{eq:ana}).  }
\label{fig:model_2}
\end{figure}

\begin{figure*}
\includegraphics[width=0.95\textwidth]{sample_model.pdf}
\caption{Sample model result with $r_c =5\times10^{17}$cm and $F=10^{10.2}$cm$^{-2}$s$^{-1}$.
The profiles are centered on the ionization front defined as where the ionization fraction is
50\%.  Note that we do not see a rise in the neutral HI fraction right
after the IF because this is where the dynamical model stops and the \textsc{Cloudy} model
begins. The \textsc{Cloudy} model starts with the HI fraction being $1$.  }
\label{fig:model_sample}
\end{figure*}

\subsection{Synthetic Br$\gamma$ and H$_2$ Spatial Line Profiles}
{Synthetic observations play a key role in interpreting real observations. For a recent review of applications of synthetic observations, see \citet{haworth}. In this investigation, we generate synthetic spatial emission profiles from the
photoevaporative flow models described above.} We create the profile for each line separately
and then stitch them together. To create the Br$\gamma$ profile, we simply square the HII
density profile calculated in the model. Because we are not interested in the actual Br$\gamma$ emissivity, we
do not need to worry about the recombination coefficient. In doing this, we ignore
the dependence of the recombination coefficient on the temperature but that is a minor
effect compared to the variations in HII density.

To create the H$_2$ profile, we use
the photoionization code \textsc{Cloudy}, version 13.03, described by \citet{ferland2013}. \textsc{Cloudy} works
by providing it a density law and input spectrum and it does the necessary radiative
transfer to predict the H$_{2}$ emissivity. \citet{bertoldi1996} found that, after a transitory initial phase,
the dissociation front (DF) is trapped within the shocked layer. This means that the PDR,
where the HI exists and where the 2.12$\mu$m emission is coming from, is situated in the
shocked layer. Hence, we model the PDR in \textsc{Cloudy} by assuming a constant density of
$n_{\rm I}$ found from the dynamical models. The shape of the input spectrum is that
of an O6 star from the TLUSTY stellar atmosphere models \citep{hubeny1995} with all of the
Lyman continuum photons removed. This is only an approximation as the photoevaporative flow will
modify the spectrum more than just extinguishing the Lyc photons, however it should be enough
to understand general trends in the offsets.

The models are also approximate in the sense
that \textsc{Cloudy} does static modeling but the PDR is actually evolving as hydrogen advects
across it at the speed of the dissociation front. \citet{bertoldi1996} find that non-equilibrium effects in
the PDR, including transient photochemistry, are important when the incident flux and the density of
the cloud is low. From their Figure~4, the regions we consider that have strong incident
flux like that characteristic of Carina and Cygnus should be well-matched with equilibrium PDR models (like \textsc{Cloudy}) though objects with
weaker incident flux like in IC 1396 might not be. Despite this limitation, we use \textsc{Cloudy} to get a rough
understanding of the general trends in the offsets.
The models run deep into the molecular region of the clouds, and we adopt \textsc{Cloudy}'s default
chemical abundances for an HII region, which are average values of the
abundances in the Orion Nebula \citet{rubin1991,baldwin1991,osterbrock1992,rubin1993}, with
a metallicity of 0.4 solar. The models also include dust grains with
the size distribution and abundances of grains found in the Orion Nebula. Additionally, the models
include a cosmic ray flux characteristic of the galactic background taken from \citet{indrilio2007}. We combine
the two profiles together simply by stitching the H$_2$ profile onto the location
where the dynamical models begin (with $F_0 = 10^{-8}F_{\rm inc}$ and $x_0 = 10^{-8}$) . Both the Br$\gamma$ and
H$_2$ profiles are then blurred to the same level of seeing as the NIR
images by convolving them with the appropriate Gaussian. A sample model result is shown in
Figure \ref{fig:model_sample}. The density profile shows the $1/r^{2}$ fall-of in the photoionized flow,
the large density contrast across the IF, and the constant density through the PDR. The
offset is then easily measured as the distance between the peak Br$\gamma$ and H$_2$ emission.

{Finally, it is worth noting that the models are 1D whereas the real objects are 3D and so these spatial profiles will not capture the effects of projecting the 3D emission profiles along the line of sight. This will have the effect of blurring out the profiles since sight lines will probe ionized material at different locations in the photoevaporative flow and neutral material. However, since the radii of curvature of the objects we consider are generally $\gtrsim 5$ times the offset size, the effect of projection along the line of sight will be relatively small on the measured offset size. Accurately projecting along the line of sight would require a radiative transfer calculation of the H$_2$ through the pillar which is beyond the scope of the current project. For other successful examples of using 1D models to interpret observed spatial line profiles in photo-evaporative flows, see, e.g., \citet{sankrit, mcleod2016}. However, 1D models might lead to erroneous interpretations of spatially resolved line ratios as discussed in \citet{ercolano2011}.}

Figure \ref{fig:offset_model} shows the results of these models for a variety of values of $r_c$ and $F$. The model offsets always lie within roughly a factor of two of
$10^{17}$cm even over a wide range in parameters, in agreement with the observations.
With the exception of the left edge,
the offsets decrease with increasing $F$ and decreasing $r_c$. This trend arises from
the increase in density in the shocked layer that occurs by increasing $F$ or
decreasing $r_c$ (see Figure \ref{fig:model_2}). Increasing the density leads to a thinner PDR
because the thickness of the PDR is set by the column density of hydrogen and
a thinner PDR leads to a smaller offset. The decrease in offsets for log(F) $\lesssim$ 10.4
in Figure \ref{fig:offset_model} is caused by the incident radiation being so weak that a distinct
layer of HI is not formed in the \textsc{Cloudy} models and so the H$_2$ emission
comes from very close to the IF and the offset declines. At these low
incident fluxes, the static models of the PDR are probably inadequate as the advection of
hydrogen through the PDR will start to matter \citep{bertoldi1996}.  It is possible for the ionization
and dissociation fronts to merge under certain conditions \citep[e.g.][]{storzer1998, henney2007},
however our static models are certainly inadequate to describe this situation.

\begin{figure}
\includegraphics[width=0.47\textwidth]{offset_theory.pdf}
\caption{Predicted model offsets between Br$\gamma$ and H$_2$
for different radii of curvature and incident fluxes. }
\label{fig:offset_model}
\end{figure}

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

In this work, we present sky-subtracted large-scale mosaics in H$_{2}$ and Br$\gamma$
for three regions of massive star formation: Cyg~OB2, Carina, and IC~1396. The three
regions span roughly three orders of magnitude in ionizing flux and so provide
a powerful means to investigate the physics that underlies the
spatial profiles and offsets we observe in the dozens of PDR fronts that populate these regions.

Each of the three star forming regions contains PDRs that have
a wide variety of geometrical shapes,
ranging from more or less spherical globules to highly-elongated cylindrical pillars.
To understand how the shapes of the interfaces affect what is observed in the images, we
fit power laws to the spatial profiles of continuum-corrected Br$\gamma$ flux along
several hundred lines of sight and grouped the results according to the overall shapes of the
irradiated interfaces. The results confirm model predictions
that the shape of the object strongly influences the spatial profile of the emission in a
photoevaporative flow. Br$\gamma$ fluxes from cylindrically-shaped
clouds have a relatively flat radial power law index of $-2.54\pm0.16$,
while spherically-shaped clouds decline more rapidly with distance, with
an average power-law index of $-3.45\pm0.22$. The globule power-law exponents agree well with the expected values between $-3$ and $-4$, although the pillar power-law exponents are somewhat higher than the expected values between $-1$ and $-2$, predicted for photoevaporative flows of constant
velocity. Broadly speaking, the results are consistent with model expectations that more spherical
clouds should exhibit sharper radial decline in the Br$\gamma$ emission
as compared with cylindrical clouds.

We observe spatial offsets
between the peaks of the H$_{2}$ and Br$\gamma$ emission throughout
each of the three star formation regions, always in the sense that
the H$_2$ emission is located deeper into the molecular cloud, underneath
a layer of Br$\gamma$ emission that traces a photoevaporative flow.
By extracting several hundreds of
slices across a wide variety of distinct PDR interfaces, we
found that the spatial offsets between these two lines are
always on the order of $10^{17}$~cm
for all the objects in all three regions.

After estimating the incident ionizing flux, $F$,
received by each object and measuring an
approximate radius of curvature $r_c$ for each object,
we investigated how these two variables correlate with the
observed spatial offsets between Br$\gamma$ and H$_2$.
There is a significant increase in the spatial offset as
$r_c$ increases. A weak trend of decreasing offset with $F$
is also observed.
Individual geometrical complications introduce significant scatter to these trends.

To investigate the physical processes that determine why the offsets form and what
sets their sizes, a simple theoretical model of Br$\gamma$ emission from a
photoevaporative flow is developed by solving the time-independent Euler and ionization equations.
These density profiles are then used in the numerical code \textsc{Cloudy}
to predict synthetic spatial line profiles for the H$_2$ emission.
The models reproduce all the main results in the data, including the size
of the spatial offsets between H$_2$ and Br$\gamma$, the
positive correlations of the offsets with the radius of curvature of the object,
and even the weak negative correlation of the offset with the amount of incident ionizing
flux.

Overall, these observations support the general picture of molecular
clouds being photoevaporated by young massive stars and lend some confidence that the basic physics that underlies these systems is understood. However, it is clear that
the geometries of the interfaces can be very complex, and this complexity has a significant
effect upon what is observed. Both the H$_2$ and Br$\gamma$ lines will be easily observable
at high spatial resolution in PDRs with the James Webb Space Telescope, and all of the
analyses done here will be possible to extend to much smaller spatial scales as long as
continuum images are also acquired to enable accurate background subtraction of the
ambient reflected starlight.

\acknowledgments
This research was supported by Department of Energy's
NNSA division under the NLUF grant DE-FOA-0001109,
and by NASA's National Radio Astronomy student support grant SOSPA3-018.
Financial support for this project was provided by the Department of Energy
grant GR523126.
Based on observations at Cerro Tololo Inter-American Observatory, National Optical Astronomy Observatory (NOAO 11A-0234; PI: P. Hartigan), which is operated by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the National Science Foundation.
Based on observations at Kitt Peak National Observatory, National Optical Astronomy Observatory (NOAO 12B-0166; PI: P. Hartigan), which is operated by the Association of Universities for Research in Astronomy (AURA) under cooperative agreement with the National Science Foundation.

{\software{Cloudy \citep{ferland2013}, NEWFIRM Pipeline \citep{swaters2009}}}

\bibliographystyle{aasjournal}