Scale size of cometary bow shocks

N. J. T. Edberg, A. I. Eriksson, E. Vigren, H. Nilsson, H. Gunell, C. Götz, I.

Richter, P. Henri, J. de Keyser

► To cite this version:

N. J. T. Edberg, A. I. Eriksson, E. Vigren, H. Nilsson, H. Gunell, et al. Scale size of cometary bow shocks. Astronomy and Astrophysics - A&A, 2024, 682, 10.1051/0004-6361/202346566. insu-04462075

HAL Id: insu-04462075

https://insu.hal.science/insu-04462075

Submitted on 16 Feb 2024

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_1.png" /></figure>

HAL
open science

Scale size of cometary bow shocks

N. J. T. Edberg, A. I. Eriksson, E. Vigren, H. Nilsson, H. Gunell, C. Götz, I.

Richter, P. Henri, J. de Keyser

► To cite this version:

N. J. T. Edberg, A. I. Eriksson, E. Vigren, H. Nilsson, H. Gunell, et al. Scale size of cometary bow shocks. Astronomy & Astrophysics - A&A, 2024, 682, 10.1051/0004-6361/202346566. insu-04462075

HAL Id: insu-04462075

https://insu.hal.science/insu-04462075v1

Submitted on 16 Feb 2024

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Scale size of cometary bow shocks

N. J. T. Edberg¹, A.I. Eriksson¹, E. Vigren¹, H. Nilsson², H. Gunell³, C. Götz⁴,
I. Richter⁵, P. Henri⁶,⁷, and J. De Keyser⁸

¹ Swedish Institute of Space Physics (IRF), Uppsala, Sweden
e-mail: ne@irfu.se

² Swedish Institute of Space Physics (IRF), Kiruna, Sweden

³ Department of Physics, Umeå University, Umeå, Sweden

⁴ Department of Mathematics, Physics and Electrical Engineering, Northumbria University, Newcastle-upon-Tyne, UK

⁵ Institut für Geophysik und extraterrestrische Physik, Technische Universität Braunschweig, Braunschweig, Germany

⁶ Laboratoire de Physique et Chimie de l'Environnement et de l'Espace, CNRS, 45071 Orléans, France

⁷ Laboratoire Lagrange, OCA, CNRS, UCA, 06304 Nice, France

⁸ Royal Belgian Institute for Space Aeronomy, BIRA-IASB, Brussels, Belgium

Received 31 March 2023 / Accepted 21 October 2023

## ABSTRACT

**Context.** In past decades, several spacecraft have visited comets to investigate their plasma environments. In the coming years, Comet Interceptor will make yet another attempt. This time, the target comet and its outgassing activity are unknown and may not be known before the spacecraft has been launched into its parking orbit, where it will await a possible interception. If the approximate outgassing rate can be estimated remotely when a target has been identified, it is desirable to also be able to estimate the scale size of the plasma environment, defined here as the region bound by the bow shock.

**Aims.** This study aims to combine previous measurements and simulations of cometary bow shock locations to gain a better understanding of how the scale size of cometary plasma environments varies. We compare these data with models of the bow shock size, and we furthermore provide an outgassing rate-dependent shape model of the bow shock. We then use this to predict a range of times and cometocentric distances for the crossing of the bow shock by Comet Interceptor, together with expected plasma density measurements along the spacecraft track.

**Methods.** We used data of the location of cometary bow shocks from previous spacecraft missions, together with simulation results from previously published studies. We compared these results with an existing model of the bow shock stand-off distance and expand on this to provide a shape model of cometary bow shocks. The model in particular includes the cometary outgassing rate, but also upstream solar wind conditions, ionisation rates, and the neutral flow velocity.

**Results.** The agreement between the gas-dynamic model and the data and simulation results is good in terms of the stand-off distance of the bow shock as a function of the outgassing rate. For outgassing rates in the range of $10^{27}-10^{31}$ s$^{-1}$, the scale size of cometary bow shocks can vary by four orders of magnitude, from about $10^2$ km to $10^6$ km, for an ionisation rate, flow velocity, and upstream solar wind conditions typical of those at 1 AU. The proposed bow shock shape model shows that a comet plasma environment can range in scale size from the plasma environment of Mars to about half of that of Saturn.

**Conclusions.** The model-data agreement allows for the planning of upcoming spacecraft comet encounters, such as that of Comet Interceptor, when a target has been identified and its outgassing rate is determined. We conclude that the time a spacecraft can spend within the plasma environment during a flyby can range from minutes to days, depending on the comet that is visited and on the flyby speed. However, to capture most of the comet plasma environment, including pick-up ions and upstream plasma waves, and to ensure the highest possible scientific return, measurements should still start well upstream of the expected bow shock location. From the plasma perspective, the selected target should preferably be an active comet with the lowest possible flyby velocity.

**Key words.** comets: general – plasmas

## 1. Introduction

The plasma environments of active comets are some of the largest structures in the solar system. However, studying them in situ and obtaining successful measurements requires meticulously planned spacecraft missions optimised for capturing often fleeting opportunities. Hence, it is crucial to determine in advance the spatial extent of the region that needs to be studied. This paper focuses on investigating the scale size of the outermost plasma region of comets, delineated by the bow shock, and its dependence on the cometary outgassing rate. Although plasma waves and cometary ions may still exist beyond this boundary, we define the bow shock as the outermost clearly identifiable plasma boundary of the cometary

plasma environment. Understanding this boundary provides a first-order approximation for designing the operational timeline of future cometary flyby missions, such as the upcoming Comet Interceptor of the European Space Agency (Snodgrass & Jones 2019; Jones et al. 2022).

Cometary flybys and visits have been conducted several times since the initial endeavour in 1985 when the International Cometary Explorer (ICE) flew past comet 21P/Giacobini-Zinner (Rosenvinge et al. 1986; Bame et al. 1986). In 1986, multiple spacecraft from Europe, Japan, and the Soviet Union, collectively known as the Halley Armada, visited comet 1P/Halley: Giotto (Reinhard 1986), Sakigake and Suisei (Hirao & Itoh 1986), and Vega 1 and 2 (Sagdeev et al. 1986). Although the ICE spacecraft also flew past 1P/Halley, along with
# Scale size of cometary bow shocks

N. J. T. Edberg¹, A.I. Eriksson¹, E. Vigren¹, H. Nilsson², H. Gunell³, C. Götz⁴,
I. Richter⁵, P. Henri⁶,⁷, and J. De Keyser⁸

¹ Swedish Institute of Space Physics (IRF), Uppsala, Sweden
e-mail: ne@irfu.se

² Swedish Institute of Space Physics (IRF), Kiruna, Sweden

³ Department of Physics, Umeå University, Umeå, Sweden

⁴ Department of Mathematics, Physics and Electrical Engineering, Northumbria University, Newcastle-upon-Tyne, UK

⁵ Institut für Geophysik und extraterrestrische Physik, Technische Universität Braunschweig, Braunschweig, Germany

⁶ Laboratoire de Physique et Chimie de l'Environnement et de l'Espace, CNRS, 45071 Orléans, France

⁷ Laboratoire Lagrange, OCA, CNRS, UCA, 06304 Nice, France

⁸ Royal Belgian Institute for Space Aeronomy, BIRA-IASB, Brussels, Belgium

Received 31 March 2023 / Accepted 21 October 2023

## ABSTRACT

**Context.** In past decades, several spacecraft have visited comets to investigate their plasma environments. In the coming years, Comet Interceptor will make yet another attempt. This time, the target comet and its outgassing activity are unknown and may not be known before the spacecraft has been launched into its parking orbit, where it will await a possible interception. If the approximate outgassing rate can be estimated remotely when a target has been identified, it is desirable to also be able to estimate the scale size of the plasma environment, defined here as the region bound by the bow shock.

**Aims.** This study aims to combine previous measurements and simulations of cometary bow shock locations to gain a better understanding of how the scale size of cometary plasma environments varies. We compare these data with models of the bow shock size, and we furthermore provide an outgassing rate-dependent shape model of the bow shock. We then use this to predict a range of times and cometocentric distances for the crossing of the bow shock by Comet Interceptor, together with expected plasma density measurements along the spacecraft track.

**Methods.** We used data of the location of cometary bow shocks from previous spacecraft missions, together with simulation results from previously published studies. We compared these results with an existing model of the bow shock stand-off distance and expand on this to provide a shape model of cometary bow shocks. The model in particular includes the cometary outgassing rate, but also upstream solar wind conditions, ionisation rates, and the neutral flow velocity.

**Results.** The agreement between the gas-dynamic model and the data and simulation results is good in terms of the stand-off distance of the bow shock as a function of the outgassing rate. For outgassing rates in the range of $10^{27}-10^{31}$ s$^{-1}$, the scale size of cometary bow shocks can vary by four orders of magnitude, from about $10^2$ km to $10^6$ km, for an ionisation rate, flow velocity, and upstream solar wind conditions typical of those at 1 AU. The proposed bow shock shape model shows that a comet plasma environment can range in scale size from the plasma environment of Mars to about half of that of Saturn.

**Conclusions.** The model-data agreement allows for the planning of upcoming spacecraft comet encounters, such as that of Comet Interceptor, when a target has been identified and its outgassing rate is determined. We conclude that the time a spacecraft can spend within the plasma environment during a flyby can range from minutes to days, depending on the comet that is visited and on the flyby speed. However, to capture most of the comet plasma environment, including pick-up ions and upstream plasma waves, and to ensure the highest possible scientific return, measurements should still start well upstream of the expected bow shock location. From the plasma perspective, the selected target should preferably be an active comet with the lowest possible flyby velocity.

**Key words.** comets: general – plasmas

## 1. Introduction

The plasma environments of active comets are some of the largest structures in the solar system. However, studying them in situ and obtaining successful measurements requires meticulously planned spacecraft missions optimised for capturing often fleeting opportunities. Hence, it is crucial to determine in advance the spatial extent of the region that needs to be studied. This paper focuses on investigating the scale size of the outermost plasma region of comets, delineated by the bow shock, and its dependence on the cometary outgassing rate. Although plasma waves and cometary ions may still exist beyond this boundary, we define the bow shock as the outermost clearly identifiable plasma boundary of the cometary

plasma environment. Understanding this boundary provides a first-order approximation for designing the operational timeline of future cometary flyby missions, such as the upcoming Comet Interceptor of the European Space Agency (Snodgrass & Jones 2019; Jones et al. 2022).

Cometary flybys and visits have been conducted several times since the initial endeavour in 1985 when the International Cometary Explorer (ICE) flew past comet 21P/Giacobini-Zinner (Rosenvinge et al. 1986; Bame et al. 1986). In 1986, multiple spacecraft from Europe, Japan, and the Soviet Union, collectively known as the Halley Armada, visited comet 1P/Halley: Giotto (Reinhard 1986), Sakigake and Suisei (Hirao & Itoh 1986), and Vega 1 and 2 (Sagdeev et al. 1986). Although the ICE spacecraft also flew past 1P/Halley, along with
<figure><img src="image_2.png" /></figure>

**Fig. 1.** Time series of magnetic field data during previous cometary flybys. The vertical dashed lines indicate the times of bow shock crossings, and the dotted lines represent the time of closest approach. The flybys (panels) are ordered in decreasing time spent within the respective plasma environment. Vega 2 did not observe the actual bow shock crossing with its magnetometer, and the distance of Sakigake at the flyby was too large for it to enter the plasma environment. The Suisei spacecraft did not carry a magnetometer and is therefore not included in this plot.

Sakigake, it remained too far from the coma to enter it. After the 1P/Halley flyby, Giotto successfully encountered comet 26P/Grigg-Skjellerup in 1992 (e.g., Johnstone et al. 1993). A decade later, in 2001, Deep Space 1 performed a flyby of comet 19P/Borrelly (Young et al. 2004). Another decade later, in 2014, the Rosetta spacecraft entered the coma of comet 67P/Churyumov-Gerasimenko (Glassmeier et al. 2007; Goetz et al. 2022) and conducted observations for 2 yr. Although a few other spacecraft have made close-up measurements at comets, none of them carried instruments specifically designed for studying the plasma environment.

## 2. Observations of cometary bow shocks

The variety of comet observations is summarised in Fig. 1, which presents time series of magnetic field strength measurements from seven previous cometary flybys. The top four panels display data from comet 1P/Halley, followed by data from comets 19P, 21P, and 26P. Each dataset is centred around the time of closest approach for each flyby and covers a time span of ±6 h, including the inbound and outbound bow shock crossings. It is worth noting that there has been some discussion in the literature regarding whether all these bow shock crossings can indeed be classified

as shocks. In some cases, they are referred to as “bow waves” due to their wide width, gradual changes in plasma parameters, and a Mach number below one, as seen in the inbound crossing at 26P (Coates et al. 1997). However, for simplicity, we adopt the term “bow shock” throughout this paper. The panels are arranged in descending order based on the time spent within the bow shock, ranging from approximately 7.5 h for Giotto at 1P/Halley to around 1 h for Giotto at comet 26P. The cadence of the magnetic field measurements varies across the panels, and the flyby velocity differs in each case.

The trajectories of some of the previous spacecraft that have flown past comets and observed their respective bow shocks are depicted in Fig. 2. These trajectories are illustrated in the cometary solar ecliptic (CSEQ) reference frame, which employs a right-handed system. In this frame, the x-axis extends from the comet nucleus towards the Sun, the z-axis repre-sents the component of the solar north pole of date orthogonal to the x-axis, and the y-axis completes the reference frame. Each reference frame is centred on its respective comet. It is worth noting that the flyby trajectories of spacecraft visiting comets have generally followed a plane close to the termina-tor plane, allowing for the sampling of cometary bow shocks in this region.
<figure><img src="image_3.png" /></figure>

**Fig. 1.** Time series of magnetic field data during previous cometary flybys. The vertical dashed lines indicate the times of bow shock crossings, and the dotted lines represent the time of closest approach. The flybys (panels) are ordered in decreasing time spent within the respective plasma environment. Vega 2 did not observe the actual bow shock crossing with its magnetometer, and the distance of Sakigake at the flyby was too large for it to enter the plasma environment. The Suisei spacecraft did not carry a magnetometer and is therefore not included in this plot.

Sakigake, it remained too far from the coma to enter it. After the 1P/Halley flyby, Giotto successfully encountered comet 26P/Grigg-Skjellerup in 1992 (e.g., Johnstone et al. 1993). A decade later, in 2001, Deep Space 1 performed a flyby of comet 19P/Borrelly (Young et al. 2004). Another decade later, in 2014, the Rosetta spacecraft entered the coma of comet 67P/Churyumov-Gerasimenko (Glassmeier et al. 2007; Goetz et al. 2022) and conducted observations for 2 yr. Although a few other spacecraft have made close-up measurements at comets, none of them carried instruments specifically designed for studying the plasma environment.

## 2. Observations of cometary bow shocks

The variety of comet observations is summarised in Fig. 1, which presents time series of magnetic field strength measurements from seven previous cometary flybys. The top four panels display data from comet 1P/Halley, followed by data from comets 19P, 21P, and 26P. Each dataset is centred around the time of closest approach for each flyby and covers a time span of ±6 h, including the inbound and outbound bow shock crossings. It is worth noting that there has been some discussion in the literature regarding whether all these bow shock crossings can indeed be classified

as shocks. In some cases, they are referred to as “bow waves” due to their wide width, gradual changes in plasma parameters, and a Mach number below one, as seen in the inbound crossing at 26P (Coates et al. 1997). However, for simplicity, we adopt the term “bow shock” throughout this paper. The panels are arranged in descending order based on the time spent within the bow shock, ranging from approximately 7.5 h for Giotto at 1P/Halley to around 1 h for Giotto at comet 26P. The cadence of the magnetic field measurements varies across the panels, and the flyby velocity differs in each case.

The trajectories of some of the previous spacecraft that have flown past comets and observed their respective bow shocks are depicted in Fig. 2. These trajectories are illustrated in the cometary solar ecliptic (CSEQ) reference frame, which employs a right-handed system. In this frame, the x-axis extends from the comet nucleus towards the Sun, the z-axis repre-sents the component of the solar north pole of date orthogonal to the x-axis, and the y-axis completes the reference frame. Each reference frame is centred on its respective comet. It is worth noting that the flyby trajectories of spacecraft visiting comets have generally followed a plane close to the termina-tor plane, allowing for the sampling of cometary bow shocks in this region.
<figure><img src="image_3.png" /><figcaption>Fig. 2. Some previous cometary flybys shown in cylindrical coordinates in the CSEQ reference frame for each respective comet. The bow shock models are shown with solid black lines and are produced as described in the text. The locations of the spacecraft bow shock crossings are indicated by symbols along the tracks. We took the average between the inbound and the outbound locations of the shock crossings to fit to the model.</figcaption></figure>

In most cases, the bow shock crossings are asymmetric around the closest approach. This is due either to a flyby trajectory that is tilted with respect to the comet-Sun line and/or due to an asymmetric plasma environment arising from, for instance, the direction of the interplanetary magnetic field, solar wind flow shears, finite ion gyroradius effects, or a spatially nonuniform outgassing. The solar wind conditions can also have changed

from the time of the inbound to the outbound leg, resulting in a displacement of the bow shock location, as shown in Behar et al. (2022); Behar & Henri (2023). The asymmetry can also be caused by aberration of the solar wind direction by an angle determined by the cometary velocity vector and the mean flow direction of the solar wind: in the comet frame of reference, the solar wind would appear to come from a direction different to that of the Sun, and in the CSEQ reference frame, the bow shock is tilted in the direction opposite of the comet movement.

The scale size of the encountered plasma environments, defined here as the region bound by the bow shock, ranges from about 10⁶ km in the case of 1P/Halley to a few times 10⁴ km in the cases of comet 19P/Borrelly, 21P/Giacobini-Zinner, and 26P/Grigg-Skjellerup (note the change in scale from the first panel in Fig. 2 to the following three panels). The location in space of the bow shock crossings, indicated by the different symbols in each panel, are determined from primarily magnetic field measurements, and we used the times and locations already determined by the numerous previous studies of these crossings (e.g. Galeev 1985; Gringauz et al. 1986; Mukai et al. 1986; Jones et al. 1986; Johnstone et al. 1993; Richter et al. 2011; Mazelle et al. 2004). Some times differ slightly between studies, depending on interpretation and the data set being used, but for our purpose here, this discrepancy is insignificant. The shape of the bow shocks shown in Fig. 2 was not included in the previous studies, but was added here, as we describe next.

The crossings are mostly in the terminator plane, and little is therefore known of the subsolar bow shock location and consequently also of the overall shape of the bow shocks, at least from measurements. Flammer & Mendis (1993) used the concept of flaring to extrapolate the terminator crossings of the bow shock of 26P/Grigg-Skjellerup to the subsolar point, that is, they assumed that the bow shock shape was a paraboloid and that the terminator distance was, in that particular case, 2.15 times the subsolar distance.

At many planets, both magnetised and un-magnetised, the shapes of the bow shocks have been found empirically from many months or years of orbiting spacecraft data to be well described by conic sections of the form

$$r = \frac{L}{1 + \epsilon \cos \theta'} \quad (1)$$

where $(r, \theta)$ are the polar coordinates in an aberrated CSEQ frame $(X, Y, Z)$ such that the mean solar wind flow is in the $-X$-direction in the comet frame of reference (to account for the fact that the cometary transverse velocity $v_y$ is non-negligible compared to the solar wind radial velocity $v_r$, i.e. the CSEQ frame is rotated about the $z$-axis with an angle $-\arctan(v_y/v_r)$), $L$ is the semi-latus rectum, and $\epsilon$ is the eccentricity (e.g. Farris et al. 1991; Huddleston et al. 1998; Masters et al. 2008; Edberg et al. 2008; Martinecz et al. 2009). For Mars and Venus, which are relatively large obstacles to the solar wind compared to their bow shock sizes, the conic sections are centred at $X, Y, Z = (X_0, 0, 0)$, where $X_0$ is a free model parameter. For bow shock models at other planets that are small compared to their bow shocks, such as Earth, Saturn, and Jupiter, the conic sections are centred at $(0, 0, 0)$. The shapes shown in Fig. 2 therefore assume a conic section shape because that is generally what has been observed elsewhere in the solar system.

The conic section model has two free model parameters $L$ and $\epsilon$ ($X_0$ is assumed equal to 0 for comets as the bow shock distance is much larger than the comet nucleus), and we cannot know the proper shape of the bow shock. We therefore simply
<figure><img src="image_4.png" /><figcaption>Fig. 2. Some previous cometary flybys shown in cylindrical coordinates in the CSEQ reference frame for each respective comet. The bow shock models are shown with solid black lines and are produced as described in the text. The locations of the spacecraft bow shock crossings are indicated by symbols along the tracks. We took the average between the inbound and the outbound locations of the shock crossings to fit to the model.</figcaption></figure>

In most cases, the bow shock crossings are asymmetric around the closest approach. This is due either to a flyby trajectory that is tilted with respect to the comet-Sun line and/or due to an asymmetric plasma environment arising from, for instance, the direction of the interplanetary magnetic field, solar wind flow shears, finite ion gyroradius effects, or a spatially nonuniform outgassing. The solar wind conditions can also have changed

from the time of the inbound to the outbound leg, resulting in a displacement of the bow shock location, as shown in Behar et al. (2022); Behar & Henri (2023). The asymmetry can also be caused by aberration of the solar wind direction by an angle determined by the cometary velocity vector and the mean flow direction of the solar wind: in the comet frame of reference, the solar wind would appear to come from a direction different to that of the Sun, and in the CSEQ reference frame, the bow shock is tilted in the direction opposite of the comet movement.

The scale size of the encountered plasma environments, defined here as the region bound by the bow shock, ranges from about 10⁶ km in the case of 1P/Halley to a few times 10⁴ km in the cases of comet 19P/Borrelly, 21P/Giacobini-Zinner, and 26P/Grigg-Skjellerup (note the change in scale from the first panel in Fig. 2 to the following three panels). The location in space of the bow shock crossings, indicated by the different symbols in each panel, are determined from primarily magnetic field measurements, and we used the times and locations already determined by the numerous previous studies of these crossings (e.g. Galeev 1985; Gringauz et al. 1986; Mukai et al. 1986; Jones et al. 1986; Johnstone et al. 1993; Richter et al. 2011; Mazelle et al. 2004). Some times differ slightly between studies, depending on interpretation and the data set being used, but for our purpose here, this discrepancy is insignificant. The shape of the bow shocks shown in Fig. 2 was not included in the previous studies, but was added here, as we describe next.

The crossings are mostly in the terminator plane, and little is therefore known of the subsolar bow shock location and consequently also of the overall shape of the bow shocks, at least from measurements. Flammer & Mendis (1993) used the concept of flaring to extrapolate the terminator crossings of the bow shock of 26P/Grigg-Skjellerup to the subsolar point, that is, they assumed that the bow shock shape was a paraboloid and that the terminator distance was, in that particular case, 2.15 times the subsolar distance.

At many planets, both magnetised and un-magnetised, the shapes of the bow shocks have been found empirically from many months or years of orbiting spacecraft data to be well described by conic sections of the form

$$r = \frac{L}{1 + \epsilon \cos \theta'} \quad (1)$$

where $(r, \theta)$ are the polar coordinates in an aberrated CSEQ frame $(X, Y, Z)$ such that the mean solar wind flow is in the $-X$-direction in the comet frame of reference (to account for the fact that the cometary transverse velocity $v_y$ is non-negligible compared to the solar wind radial velocity $v_r$, i.e. the CSEQ frame is rotated about the $z$-axis with an angle $-\arctan(v_y/v_r)$), $L$ is the semi-latus rectum, and $\epsilon$ is the eccentricity (e.g. Farris et al. 1991; Huddleston et al. 1998; Masters et al. 2008; Edberg et al. 2008; Martinecz et al. 2009). For Mars and Venus, which are relatively large obstacles to the solar wind compared to their bow shock sizes, the conic sections are centred at $X, Y, Z = (X_0, 0, 0)$, where $X_0$ is a free model parameter. For bow shock models at other planets that are small compared to their bow shocks, such as Earth, Saturn, and Jupiter, the conic sections are centred at $(0, 0, 0)$. The shapes shown in Fig. 2 therefore assume a conic section shape because that is generally what has been observed elsewhere in the solar system.

The conic section model has two free model parameters $L$ and $\epsilon$ ($X_0$ is assumed equal to 0 for comets as the bow shock distance is much larger than the comet nucleus), and we cannot know the proper shape of the bow shock. We therefore simply
<table><caption>**Table 1.** Bow shock shape model parameters (see Eq. (1)) together with the bow shock terminator and subsolar stand-off distances for some planets and comets.</caption><thead><tr><th>Comet/planet</th><th>L (km)</th><th>ε</th><th>X<sub>0</sub> (km)</th><th>R<sub>BS</sub> (km)</th><th>R<sub>T</sub> (km)</th><th>Reference</th></tr></thead><tbody><tr><td>Mars</td><td>3.0 × 10<sup>3</sup></td><td>0.92</td><td>2.9 × 10<sup>3</sup></td><td>4.5 × 10<sup>3</sup></td><td>4.9 × 10<sup>3</sup></td><td>Edberg et al. (2008)</td></tr><tr><td>Venus</td><td>9.2 × 10<sup>3</sup></td><td>1.02</td><td>4.0 × 10<sup>3</sup></td><td>8.5 × 10<sup>3</sup></td><td>1.3 × 10<sup>4</sup></td><td>Martinecz et al. (2009)</td></tr><tr><td>Earth</td><td>1.6 × 10<sup>5</sup></td><td>0.81</td><td>0</td><td>8.7 × 10<sup>5</sup></td><td>1.6 × 10<sup>5</sup></td><td>Farris et al. (1991)</td></tr><tr><td>Jupiter</td><td>9.8 × 10<sup>6</sup></td><td>0.90</td><td>0</td><td>5.2 × 10<sup>6</sup></td><td>9.8 × 10<sup>6</sup></td><td>Huddleston et al. (1998)</td></tr><tr><td>Saturn</td><td>3.4 × 10<sup>6</sup></td><td>0.92</td><td>0</td><td>1.8 × 10<sup>6</sup></td><td>3.4 × 10<sup>6</sup></td><td>Masters et al. (2008)</td></tr><tr><td>1P/Halley</td><td>7.6 × 10<sup>5</sup></td><td>1.0</td><td>0</td><td>3.8 × 10<sup>5</sup></td><td>7.6 × 10<sup>5</sup></td><td>This work</td></tr><tr><td>19P/Borrelly</td><td>1.2 × 10<sup>5</sup></td><td>1.0</td><td>0</td><td>6.8 × 10<sup>4</sup></td><td>1.3 × 10<sup>5</sup></td><td>This work</td></tr><tr><td>21P/G-Z</td><td>6.8 × 10<sup>4</sup></td><td>1.0</td><td>0</td><td>3.4 × 10<sup>4</sup></td><td>6.8 × 10<sup>4</sup></td><td>This work</td></tr><tr><td>26P/G-S</td><td>2.2 × 10<sup>4</sup></td><td>1.0</td><td>0</td><td>1.1 × 10<sup>4</sup></td><td>2.2 × 10<sup>4</sup></td><td>This work</td></tr></tbody></table>

assumed that the eccentricity is equal to 1. This introduces some error, but serves our purposes of obtaining an approximate model shape that is usable in conditions that vary by several orders of magnitude. A change in the eccentricity $\epsilon$ from 1.0 to 0.9 changes the subsolar standoff distance by no more than about 10% over the entire dayside and terminator region ($\theta < 110^{\circ}$). The bow shock model parameters for the four comets in Fig. 2 are shown in Table 1, and for comparison, the model parameters for five planets are also included. The subsolar standoff distance $R_{BS}$ and the bow shock distances at the terminator $R_{T}$ are also shown.

### 3. Simulations of cometary bow shocks

In addition to measurements, a number of numerical simulations of the comet-solar wind interaction have been performed over the past decades. These simulations have been targeting different comets of varying outgassing rates for specific purposes in each case. However, when combined, they constitute a rich source of information on scale sizes of comets, which we use in a summarising way in this paper. While we refer to the individual simulation papers for the details of each study, we briefly present some key features for each simulation in the following sections.

#### 3.1. High outgassing rate (>10<sup>29</sup> s<sup>-1</sup>)

As the scale sizes of all comets with this high outgassing rate are much larger than the ion gyroradii, magnetohydrodynamical (MHD) codes have exclusively been used to simulate their plasma environment. These codes are also more computationally efficient than hybrid or particle-in-cell simulations. Comet 1P/Halley had an outgassing rate of approximately $5 \times 10^{29}$ s<sup>-1</sup> when the armada of spacecraft flew past it in March 1986 at a heliocentric distance of about 0.89 AU. Subsequent simulations of the solar wind interaction with comet 1P/Halley were carried out by Benna & Mahaffy (2007) and Rubin et al. (2009). While no other comet with such a strong outgassing rate magnitude has been visited at around 1 AU, other comets closer to the Sun have been observed remotely, and their outgassing rate was estimated to be of the same order of magnitude as that of 1P/Halley. The non-periodic Comet C/2006 P1/McNaught had its most recent perihelion at 0.17 AU in 2007. Shou et al. (2015) simulated the plasma environment of this comet as it approached the Sun. They simulated eight different outgassing rates in the range of $7.5 \times 10^{28}–5.5 \times 10^{31}$ s<sup>-1</sup>, corresponding to heliocentric distances from 1.7 AU to 0.17 AU. Concerning the solar wind interaction with the Sun-grazing comets, another study was carried out for comet C/2011 W3 Lovejoy when that comet was at a

mere distance of $6 \times 10^{-3}$ AU from the Sun, and the outgassing rate was estimated to be $2.6 \times 10^{30}$ s<sup>-1</sup> (Jia et al. 2014).

#### 3.2. Medium outgassing rate (10<sup>28</sup>–10<sup>29</sup> s<sup>-1</sup>)

As the scale sizes of comets in this outgassing regime drop to approach the scale of the ion gyroradius, simulations are performed using either MHD codes or hybrid codes. Comets 19P/Borrelly and 21P/Giacobini-Zinner both had an outgassing rate of approximately $10^{28}$ s<sup>-1</sup> and were visited by the Deep Space 1 and the ICE spacecraft, respectively. The plasma environment of 19P/Borrelly was simulated by Delamere (2006) and Jia et al. (2008) using a hybrid code and a MHD code, respectively. The plasma environment of 19P was observed in spacecraft data to be highly asymmetric, with the bow shock location much farther out on the inbound (northern) leg compared to the outbound (southern) leg. This was initially suggested to be due to an asymmetric outgassing rate with a higher outgassing rate on the northern side of the nucleus (Young et al. 2004; Schleicher et al. 2003). However, Delamere (2006) could show that this asymmetry could arise due to the pick-up ion gyromotion depending on magnetic field direction and magnitude. To the best of our knowledge, there has not been any simulation of comet 21P, which would otherwise have fitted in this section.

Furthermore, Comet 2P/Encke had an estimated outgassing rate in the same range and was simulated by Jia et al. (2009) using the same code, but now in the context of it being impacted by a coronal mass ejection (CME) that caused a tail-disconnection event. The interaction with nominal solar wind conditions was also included in that study, and we use the values from that run here.

#### 3.3. Low outgassing rate (<10<sup>28</sup> s<sup>-1</sup>)

In this outgassing regime, it is difficult for MHD codes to obtain accurate results because kinetic effects in the plasma environment become stronger, and hybrid codes or particle-in-cell codes are more appropriate. Comet 26P/Grigg-Skjellerup, the second target of Giotto, was a low-activity comet with $Q < 10^{28}$ s<sup>-1</sup> at 1 AU, and similar to comet 67P/Churyumov-Gerasimenko when 67P was close to perihelion. As the Rosetta spacecraft escorted comet 67P from ~3.5 AU to perihelion at ~1.2 AU, the outgassing ranged from about $10^{26}$ s<sup>-1</sup> to $5 \times 10^{28}$ s<sup>-1</sup>, that is, it changed from weak to medium activity. Benna & Mahaffy (2006) used an MHD code, despite the scale size, to simulate the plasma environment of 26P and made some early predictions for what Rosetta might encounter at 67P.
<table><caption>**Table 1.** Bow shock shape model parameters (see Eq. (1)) together with the bow shock terminator and subsolar stand-off distances for some planets and comets.</caption><thead><tr><th>Comet/planet</th><th>L (km)</th><th>ε</th><th>X<sub>0</sub> (km)</th><th>R<sub>BS</sub> (km)</th><th>R<sub>T</sub> (km)</th><th>Reference</th></tr></thead><tbody><tr><td>Mars</td><td>3.0 × 10<sup>3</sup></td><td>0.92</td><td>2.9 × 10<sup>3</sup></td><td>4.5 × 10<sup>3</sup></td><td>4.9 × 10<sup>3</sup></td><td>Edberg et al. (2008)</td></tr><tr><td>Venus</td><td>9.2 × 10<sup>3</sup></td><td>1.02</td><td>4.0 × 10<sup>3</sup></td><td>8.5 × 10<sup>3</sup></td><td>1.3 × 10<sup>4</sup></td><td>Martinecz et al. (2009)</td></tr><tr><td>Earth</td><td>1.6 × 10<sup>5</sup></td><td>0.81</td><td>0</td><td>8.7 × 10<sup>5</sup></td><td>1.6 × 10<sup>5</sup></td><td>Farris et al. (1991)</td></tr><tr><td>Jupiter</td><td>9.8 × 10<sup>6</sup></td><td>0.90</td><td>0</td><td>5.2 × 10<sup>6</sup></td><td>9.8 × 10<sup>6</sup></td><td>Huddleston et al. (1998)</td></tr><tr><td>Saturn</td><td>3.4 × 10<sup>6</sup></td><td>0.92</td><td>0</td><td>1.8 × 10<sup>6</sup></td><td>3.4 × 10<sup>6</sup></td><td>Masters et al. (2008)</td></tr><tr><td>1P/Halley</td><td>7.6 × 10<sup>5</sup></td><td>1.0</td><td>0</td><td>3.8 × 10<sup>5</sup></td><td>7.6 × 10<sup>5</sup></td><td>This work</td></tr><tr><td>19P/Borrelly</td><td>1.2 × 10<sup>5</sup></td><td>1.0</td><td>0</td><td>6.8 × 10<sup>4</sup></td><td>1.3 × 10<sup>5</sup></td><td>This work</td></tr><tr><td>21P/G-Z</td><td>6.8 × 10<sup>4</sup></td><td>1.0</td><td>0</td><td>3.4 × 10<sup>4</sup></td><td>6.8 × 10<sup>4</sup></td><td>This work</td></tr><tr><td>26P/G-S</td><td>2.2 × 10<sup>4</sup></td><td>1.0</td><td>0</td><td>1.1 × 10<sup>4</sup></td><td>2.2 × 10<sup>4</sup></td><td>This work</td></tr></tbody></table>

assumed that the eccentricity is equal to 1. This introduces some error, but serves our purposes of obtaining an approximate model shape that is usable in conditions that vary by several orders of magnitude. A change in the eccentricity $\epsilon$ from 1.0 to 0.9 changes the subsolar standoff distance by no more than about 10% over the entire dayside and terminator region ($\theta < 110^{\circ}$). The bow shock model parameters for the four comets in Fig. 2 are shown in Table 1, and for comparison, the model parameters for five planets are also included. The subsolar standoff distance $R_{BS}$ and the bow shock distances at the terminator $R_{T}$ are also shown.

### 3. Simulations of cometary bow shocks

In addition to measurements, a number of numerical simulations of the comet-solar wind interaction have been performed over the past decades. These simulations have been targeting different comets of varying outgassing rates for specific purposes in each case. However, when combined, they constitute a rich source of information on scale sizes of comets, which we use in a summarising way in this paper. While we refer to the individual simulation papers for the details of each study, we briefly present some key features for each simulation in the following sections.

#### 3.1. High outgassing rate ($>10^{29} s^{-1}$)

As the scale sizes of all comets with this high outgassing rate are much larger than the ion gyroradii, magnetohydrodynamical (MHD) codes have exclusively been used to simulate their plasma environment. These codes are also more computationally efficient than hybrid or particle-in-cell simulations. Comet 1P/Halley had an outgassing rate of approximately $5 \times 10^{29} s^{-1}$ when the armada of spacecraft flew past it in March 1986 at a heliocentric distance of about 0.89 AU. Subsequent simulations of the solar wind interaction with comet 1P/Halley were carried out by Benna & Mahaffy (2007) and Rubin et al. (2009). While no other comet with such a strong outgassing rate magnitude has been visited at around 1 AU, other comets closer to the Sun have been observed remotely, and their outgassing rate was estimated to be of the same order of magnitude as that of 1P/Halley. The non-periodic Comet C/2006 P1/McNaught had its most recent perihelion at 0.17 AU in 2007. Shou et al. (2015) simulated the plasma environment of this comet as it approached the Sun. They simulated eight different outgassing rates in the range of $7.5 \times 10^{28} - 5.5 \times 10^{31} s^{-1}$, corresponding to heliocentric distances from 1.7 AU to 0.17 AU. Concerning the solar wind interaction with the Sun-grazing comets, another study was carried out for comet C/2011 W3 Lovejoy when that comet was at a

mere distance of $6 \times 10^{-3}$ AU from the Sun, and the outgassing rate was estimated to be $2.6 \times 10^{30} s^{-1}$ (Jia et al. 2014).

#### 3.2. Medium outgassing rate ($10^{28} - 10^{29} s^{-1}$)

As the scale sizes of comets in this outgassing regime drop to approach the scale of the ion gyroradius, simulations are performed using either MHD codes or hybrid codes. Comets 19P/Borrelly and 21P/Giacobini-Zinner both had an outgassing rate of approximately $10^{28} s^{-1}$ and were visited by the Deep Space 1 and the ICE spacecraft, respectively. The plasma environment of 19P/Borrelly was simulated by Delamere (2006) and Jia et al. (2008) using a hybrid code and a MHD code, respectively. The plasma environment of 19P was observed in spacecraft data to be highly asymmetric, with the bow shock location much farther out on the inbound (northern) leg compared to the outbound (southern) leg. This was initially suggested to be due to an asymmetric outgassing rate with a higher outgassing rate on the northern side of the nucleus (Young et al. 2004; Schleicher et al. 2003). However, Delamere (2006) could show that this asymmetry could arise due to the pick-up ion gyromotion depending on magnetic field direction and magnitude. To the best of our knowledge, there has not been any simulation of comet 21P, which would otherwise have fitted in this section.

Furthermore, Comet 2P/Encke had an estimated outgassing rate in the same range and was simulated by Jia et al. (2009) using the same code, but now in the context of it being impacted by a coronal mass ejection (CME) that caused a tail-disconnection event. The interaction with nominal solar wind conditions was also included in that study, and we use the values from that run here.

#### 3.3. Low outgassing rate ($<10^{28} s^{-1}$)

In this outgassing regime, it is difficult for MHD codes to obtain accurate results because kinetic effects in the plasma environment become stronger, and hybrid codes or particle-in-cell codes are more appropriate. Comet 26P/Grigg-Skjellerup, the second target of Giotto, was a low-activity comet with $Q < 10^{28} s^{-1}$ at 1 AU, and similar to comet 67P/Churyumov-Gerasimenko when 67P was close to perihelion. As the Rosetta spacecraft escorted comet 67P from ~3.5 AU to perihelion at ~1.2 AU, the outgassing ranged from about $10^{26} s^{-1}$ to $5 \times 10^{28} s^{-1}$, that is, it changed from weak to medium activity. Benna & Mahaffy (2006) used an MHD code, despite the scale size, to simulate the plasma environment of 26P and made some early predictions for what Rosetta might encounter at 67P.
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Subsolar stand-off distances of bow shocks from previous cometary flyby missions as well as from simulations of solar wind-comet interactions. The black symbols show spacecraft-measured distances, and the coloured symbols show simulations. The models of Biermann et al. (1967; dashed line) and Koenders et al. (2013; solid line) are also included, with nominal solar parameter values at 1 AU ($n_{sw} = 5 \text{ cm}^{-3}$ $u_{sw} = 400 \text{ km s}^{-1}$, $v = 7 \times 10^{-7} \text{ s}^{-1}$, and $u_n = 1 \text{ km s}^{-1}$). The grey shaded region indicates a change of a factor of two (up and down) in the solar wind dynamic pressure.</figcaption></figure>

When the outgassing of 67P was actually comparable to 26P, the Rosetta spacecraft was too close to the nucleus to cross the bow shock, and no fully developed bow shock was ever detected in the data. Consequently, no simulation-observation comparison could be performed, but a number of stand-alone simulations were indeed carried out for comet 67P. Hansen et al. (2007) performed MHD simulations for heliocentric distances of 2.7, 2.0, and 1.3 AU (corresponding to outgassing rates of $8 \times 10^{25}$, $8 \times 10^{26}$, and $5 \times 10^{27} \text{ s}^{-1}$, respectively) to prepare for the arrival of Rosetta and the foreseen different mission phases. Koenders et al. (2013) also made predictions for the arrival of Rosetta by performing simulations with five different outgassing rates in the more narrow range of $5 \times 10^{27} \text{ s}^{-1}$ to $1 \times 10^{28} \text{ s}^{-1}$, but with a more appropriate hybrid code. As the Rosetta mission entered the active escort phase, Koenders et al. (2015) provided significantly more detailed results from the inner coma from their hybrid code than what was previously published. Huang et al. (2016) used a four-fluid MHD simulation code, while Simon Wedlund et al. (2017); Alho et al. (2021) and Lindkvist et al. (2018) used yet two other hybrid codes, which we all included in our study. Simon Wedlund et al. (2017) also provided a more sophisticated quadratic surface fit to the bow shock location from their simulation results, rather than the more simple conic section. It should be mentioned that particle-in-cell

codes have also been used to simulate the plasma environment of comets (e.g. Deca et al. 2017, 2019), but not to provide a bow shock location.

## 4. Scale size of comets

The outcome of the simulations mentioned above in terms of the location of the bow shock are compiled in Fig. 3. The actual distance values are not always stated in the text in each respective study, so that in some cases, we read values off from figures. Figure 3 also includes the spacecraft-observed bow shock locations described in Sect. 2. The standoff distance ranges from about a few times $10^2$ to $10^6$ km for outgassing rates of about $10^{27}-10^{31} \text{ s}^{-1}$. The so-called infant bow shock observations reported by Gunell et al. (2018) and Goetz et al. (2021) from 67P are also included in the plot (grey dots and circles). These are not proper bow shock crossings, but rather signatures in the measured proton temperature that are thought to be indicative of a bow shock that begins to form upstream. The infant bow shock observations are sometimes made at even smaller distances and when outgassing rates are typically lower. Two models of the bow shock distance as a function of outgassing rate are also included and are discussed next.
<figure><img src="image_5.png" /><figcaption>**Fig. 3.** Subsolar stand-off distances of bow shocks from previous cometary flyby missions as well as from simulations of solar wind-comet interactions. The black symbols show spacecraft-measured distances, and the coloured symbols show simulations. The models of Biermann et al. (1967; dashed line) and Koenders et al. (2013; solid line) are also included, with nominal solar parameter values at 1 AU ($n_{sw} = 5 \text{ cm}^{-3}$ $u_{sw} = 400 \text{ km s}^{-1}$, $v = 7 \times 10^{-7} \text{ s}^{-1}$, and $u_n = 1 \text{ km s}^{-1}$). The grey shaded region indicates a change of a factor of two (up and down) in the solar wind dynamic pressure.</figcaption></figure>

When the outgassing of 67P was actually comparable to 26P, the Rosetta spacecraft was too close to the nucleus to cross the bow shock, and no fully developed bow shock was ever detected in the data. Consequently, no simulation-observation comparison could be performed, but a number of stand-alone simulations were indeed carried out for comet 67P. Hansen et al. (2007) performed MHD simulations for heliocentric distances of 2.7, 2.0, and 1.3 AU (corresponding to outgassing rates of $8 \times 10^{25}$, $8 \times 10^{26}$, and $5 \times 10^{27} \text{ s}^{-1}$, respectively) to prepare for the arrival of Rosetta and the foreseen different mission phases. Koenders et al. (2013) also made predictions for the arrival of Rosetta by performing simulations with five different outgassing rates in the more narrow range of $5 \times 10^{27} \text{ s}^{-1}$ to $1 \times 10^{28} \text{ s}^{-1}$, but with a more appropriate hybrid code. As the Rosetta mission entered the active escort phase, Koenders et al. (2015) provided significantly more detailed results from the inner coma from their hybrid code than what was previously published. Huang et al. (2016) used a four-fluid MHD simulation code, while Simon Wedlund et al. (2017); Alho et al. (2021) and Lindkvist et al. (2018) used yet two other hybrid codes, which we all included in our study. Simon Wedlund et al. (2017) also provided a more sophisticated quadratic surface fit to the bow shock location from their simulation results, rather than the more simple conic section. It should be mentioned that particle-in-cell

codes have also been used to simulate the plasma environment of comets (e.g. Deca et al. 2017, 2019), but not to provide a bow shock location.

## 4. Scale size of comets

The outcome of the simulations mentioned above in terms of the location of the bow shock are compiled in Fig. 3. The actual distance values are not always stated in the text in each respective study, so that in some cases, we read values off from figures. Figure 3 also includes the spacecraft-observed bow shock locations described in Sect. 2. The standoff distance ranges from about a few times $10^2$ to $10^6$ km for outgassing rates of about $10^{27}-10^{31} \text{ s}^{-1}$. The so-called infant bow shock observations reported by Gunell et al. (2018) and Goetz et al. (2021) from 67P are also included in the plot (grey dots and circles). These are not proper bow shock crossings, but rather signatures in the measured proton temperature that are thought to be indicative of a bow shock that begins to form upstream. The infant bow shock observations are sometimes made at even smaller distances and when outgassing rates are typically lower. Two models of the bow shock distance as a function of outgassing rate are also included and are discussed next.
The first model, shown by the dashed line in Fig. 3, is the gas-dynamic model from Biermann et al. (1967), which essentially describes the interaction between the solar wind and a one-dimensional inviscid gas flow. Using this model, Galeev et al. (1985) was able to give an expression for the bow shock subsolar standoff as a function of outgassing rate

$$R_{BS} = \frac{\nu Q m_i}{4\pi u_n n_{sw} m_{sw} u_{sw} ((\rho u_x)_{crit}^* - 1)}, \quad (2)$$

where $\nu$ is the total ionisation rate, $Q$ is the cometary outgassing rate, $u_n$ is the neutral velocity, $n_{sw}$ and $u_{sw}$ are the solar wind density and velocity, $m_i = 17$ amu is the mass of the cometary ions (assuming a water-dominated coma, which is likely at 1 AU, inside the freezing line), $m_{sw} = 1$ amu is the mass of the solar wind ions (protons), and the dimensionless quantity $(\rho u_x)_{crit}^* = 1.3$ is the so-called critical value of the normalised mass flux density. This value is required for the solar wind to slow down enough for a shock to form. The value of $(\rho u_x)_{crit}^*$ varies in the range ~1.3–1.56 in the literature, depending on the number of degrees of freedom chosen (Koenders et al. 2013), and we used the original value from Biermann et al. (1967).

Koenders et al. (2013) revisited this model and made two improvements in particular that are important for low and high outgassing rates. For weakly outgassing comets, they noted that the pick-up time $t_p$ becomes important, that is, the time from ionisation to when the ions with a gyrofrequency $\Omega_{ci}$ have reached a bulk flow velocity equal to that of the solar wind flow. A term $R_s = u_{sw}/\Omega_{ci}(\Omega_{ci}t_p - \sin(\Omega_{ci}t_p))$ was therefore included, which takes this time delay into account and thereby the distance an ion can travel before it is properly accelerated. It was also noted that the value of $t_p\Omega_{ci}/2\pi = 0.216$, determined using simulations, will vary with magnetic field strength, mass flux density, and distance to the Sun.

For high outgassing rates, Koenders et al. (2013) noted that there is in fact a maximum possible standoff distance of the bow shock because the ion mass source is limited by the rate of ionisation and the velocity of the neutrals. Hence, for arbitrary $Q$, they modified Eq. (2) to

$$R_{BS,K} = \frac{\nu Q m_i}{4\pi u_n n_{sw} m_{sw} u_{sw} ((\rho u_x)_{crit}^* - 1 + A)} - R_s, \quad (3)$$

where $A = \frac{\nu Q m_i}{4\pi u_n n_{sw} u_{sw} (u_n/v + R_n)}$ and $R_n$ is the radius of the comet nucleus. For weakly or intermediately outgassing comets, the term $A$ is negligible, but it becomes important for outgassing rates of the order of that of comet 1P/Halley. Equation (3) simplifies into $R_{BS,K} = u_n/v$ as $Q$ grows large and when $R_n$ is negligible. Equations (2) and (3) are plotted in Fig. 3 (dashed and solid line, respectively) with 1 AU values of the solar wind parameters, ionisation rate, and neutral flow velocity. The model by Koenders et al. (2013), Eq. (3), is also shown with an increase of a factor of two and a decrease in the solar wind flux (dark grey shaded region) to illustrate the sensitivity of the bow shock distance on upstream solar wind conditions. About 90% of all bow shock samples fall within this region (when the simulations by Shou et al. 2015, for a comet close to the Sun, are excluded), which serves as an estimate for the robustness of the model.

Another crude test of the robustness of this model is to list the cometary missions that did not observe a bow shock: Sakigake had a closest approach of Halley of 6.99 × 10⁶ km. ICE also passed by Halley, but at a much larger distance of 2.8 × 10⁷ km.

Rosetta did not see a bow shock during the approach phase, when the outgassing rate was below ~10²⁶ s⁻¹. All three of these missions would fall to the left of the model (solid line) in Fig. 3, where a bow shock would not be expected.

### 4.1. Outgassing-dependent bow shock shape model

Equation (3) gives the subsolar standoff distance of the bow shock, but not the overall shape of it. When we combine Eqs. (1) and (3) and assume that $\epsilon = 1$ and $X_0 = 0$, we obtain

$$R_{BS,K} = \frac{L}{1 + \cos(\theta = 0^{\circ})} \rightarrow \quad (4)$$

$$L(Q) = (1 + \epsilon) \left( \frac{\nu Q m_i}{4\pi u_n n_{sw} m_{sw} u_{sw} ((\rho u_x)_{crit}^* - 1 + A)} - R_s \right), \quad (5)$$

and the shape model of cometary bow shocks can be written as

$$r = \frac{L(Q)}{1 + \epsilon \cos \theta}. \quad (6)$$

Equation (6) allows us to plot the shape of cometary bow shocks as a function of $Q$, which is shown in Fig. 4 for outgassing rates in the range of 10²⁸–10³² s⁻¹. The two panels of this figure also include the bow shock models of the aforementioned comets and the planets Mars, Venus, Earth, Saturn and Jupiter for comparison.

### 4.2. Discussion

Evidently, the scale size of the plasma environments of comets (defined here as the region bound by the bow shock) can vary by orders of magnitude from one comet to the next, but also for a given comet along its orbit around the Sun. Depending primarily on the outgassing rate, but also on the solar wind conditions, the entire plasma environment can extend from a few times 10² km to 10⁶ km in the terminator and Sun direction for typical conditions at 1 AU. In the (ion) tail region, the plasma environment is much more extended. In any case, a comet bow shock scale size at 1 AU can be anything from a factor of 10 smaller than the bow shock of Mars to something just short of the bow shock of Saturn. The bow shocks of the magnetised planets are large because their large magnetospheres formed by a strong intrinsic magnetic field, so that a comparison to comets is somewhat incorrect, but still interesting: It provides insight into the scales of some of the largest structures of the solar system.

The bow shock stand-off model (Eq. (3)) fits the data from the various missions and simulations quite well and fits better than the original Biermann model at both high and low outgassing rates. An exception to this agreement is perhaps the simulation by Hansen et al. (2007) for the lowest outgassing rates. They used an MHD model, which might not be applicable to low outgassing rates when kinetic effects become increasingly important. The simulation of comet Lovejoy (Jia et al. 2014), which came very close to the Sun, gives the largest bow shock distance of all in this study. The simulations by Shou et al. (2015) give somewhat lower distances for the high outgassing cases than the model suggests (solid line). These simulations were made for the high ionisation rate expected close to the Sun, while the model was calculated for 1 AU conditions. A disagreement is therefore expected.

The minimum outgassing rate for which a bow shock forms (i.e. when $R_{BS}$ goes to zero) is about 3.5 × 10²⁷ s⁻¹ according to
The first model, shown by the dashed line in Fig. 3, is the gas-dynamic model from Biermann et al. (1967), which essentially describes the interaction between the solar wind and a one-dimensional inviscid gas flow. Using this model, Galeev et al. (1985) was able to give an expression for the bow shock subsolar standoff as a function of outgassing rate

$$R_{BS} = \frac{\nu Q m_i}{4\pi u_n n_{sw} m_{sw} u_{sw} ((\rho u_x)_{crit}^* - 1)}, \quad (2)$$

where $\nu$ is the total ionisation rate, $Q$ is the cometary outgassing rate, $u_n$ is the neutral velocity, $n_{sw}$ and $u_{sw}$ are the solar wind density and velocity, $m_i = 17$ amu is the mass of the cometary ions (assuming a water-dominated coma, which is likely at 1 AU, inside the freezing line), $m_{sw} = 1$ amu is the mass of the solar wind ions (protons), and the dimensionless quantity $(\rho u_x)_{crit}^* = 1.3$ is the so-called critical value of the normalised mass flux density. This value is required for the solar wind to slow down enough for a shock to form. The value of $(\rho u_x)_{crit}^*$ varies in the range ~1.3–1.56 in the literature, depending on the number of degrees of freedom chosen (Koenders et al. 2013), and we used the original value from Biermann et al. (1967).

Koenders et al. (2013) revisited this model and made two improvements in particular that are important for low and high outgassing rates. For weakly outgassing comets, they noted that the pick-up time $t_p$ becomes important, that is, the time from ionisation to when the ions with a gyrofrequency $\Omega_{ci}$ have reached a bulk flow velocity equal to that of the solar wind flow. A term $R_s = u_{sw}/\Omega_{ci}(\Omega_{ci}t_p - \sin(\Omega_{ci}t_p))$ was therefore included, which takes this time delay into account and thereby the distance an ion can travel before it is properly accelerated. It was also noted that the value of $t_p\Omega_{ci}/2\pi = 0.216$, determined using simulations, will vary with magnetic field strength, mass flux density, and distance to the Sun.

For high outgassing rates, Koenders et al. (2013) noted that there is in fact a maximum possible standoff distance of the bow shock because the ion mass source is limited by the rate of ionisation and the velocity of the neutrals. Hence, for arbitrary $Q$, they modified Eq. (2) to

$$R_{BS,K} = \frac{\nu Q m_i}{4\pi u_n n_{sw} m_{sw} u_{sw} ((\rho u_x)_{crit}^* - 1 + A)} - R_s, \quad (3)$$

where $A = \frac{\nu Q m_i}{4\pi u_n n_{sw} u_{sw} (u_n/v + R_n)}$ and $R_n$ is the radius of the comet nucleus. For weakly or intermediately outgassing comets, the term $A$ is negligible, but it becomes important for outgassing rates of the order of that of comet 1P/Halley. Equation (3) simplifies into $R_{BS,K} = u_n/v$ as $Q$ grows large and when $R_n$ is negligible. Equations (2) and (3) are plotted in Fig. 3 (dashed and solid line, respectively) with 1 AU values of the solar wind parameters, ionisation rate, and neutral flow velocity. The model by Koenders et al. (2013), Eq. (3), is also shown with an increase of a factor of two and a decrease in the solar wind flux (dark grey shaded region) to illustrate the sensitivity of the bow shock distance on upstream solar wind conditions. About 90% of all bow shock samples fall within this region (when the simulations by Shou et al. 2015, for a comet close to the Sun, are excluded), which serves as an estimate for the robustness of the model.

Another crude test of the robustness of this model is to list the cometary missions that did not observe a bow shock: Sakigake had a closest approach of Halley of 6.99 × 10⁶ km. ICE also passed by Halley, but at a much larger distance of 2.8 × 10⁷ km.

Rosetta did not see a bow shock during the approach phase, when the outgassing rate was below ~10²⁶ s⁻¹. All three of these missions would fall to the left of the model (solid line) in Fig. 3, where a bow shock would not be expected.

### 4.1. Outgassing-dependent bow shock shape model

Equation (3) gives the subsolar standoff distance of the bow shock, but not the overall shape of it. When we combine Eqs. (1) and (3) and assume that $\epsilon = 1$ and $X_0 = 0$, we obtain

$$R_{BS,K} = \frac{L}{1 + \cos(\theta = 0^{\circ})} \rightarrow \quad (4)$$

$$L(Q) = (1 + \epsilon) \left( \frac{\nu Q m_i}{4\pi u_n n_{sw} m_{sw} u_{sw} ((\rho u_x)_{crit}^* - 1 + A)} - R_s \right), \quad (5)$$

and the shape model of cometary bow shocks can be written as

$$r = \frac{L(Q)}{1 + \epsilon \cos \theta}. \quad (6)$$

Equation (6) allows us to plot the shape of cometary bow shocks as a function of $Q$, which is shown in Fig. 4 for outgassing rates in the range of 10²⁸–10³² s⁻¹. The two panels of this figure also include the bow shock models of the aforementioned comets and the planets Mars, Venus, Earth, Saturn and Jupiter for comparison.

### 4.2. Discussion

Evidently, the scale size of the plasma environments of comets (defined here as the region bound by the bow shock) can vary by orders of magnitude from one comet to the next, but also for a given comet along its orbit around the Sun. Depending primarily on the outgassing rate, but also on the solar wind conditions, the entire plasma environment can extend from a few times 10² km to 10⁶ km in the terminator and Sun direction for typical conditions at 1 AU. In the (ion) tail region, the plasma environment is much more extended. In any case, a comet bow shock scale size at 1 AU can be anything from a factor of 10 smaller than the bow shock of Mars to something just short of the bow shock of Saturn. The bow shocks of the magnetised planets are large because their large magnetospheres formed by a strong intrinsic magnetic field, so that a comparison to comets is somewhat incorrect, but still interesting: It provides insight into the scales of some of the largest structures of the solar system.

The bow shock stand-off model (Eq. (3)) fits the data from the various missions and simulations quite well and fits better than the original Biermann model at both high and low outgassing rates. An exception to this agreement is perhaps the simulation by Hansen et al. (2007) for the lowest outgassing rates. They used an MHD model, which might not be applicable to low outgassing rates when kinetic effects become increasingly important. The simulation of comet Lovejoy (Jia et al. 2014), which came very close to the Sun, gives the largest bow shock distance of all in this study. The simulations by Shou et al. (2015) give somewhat lower distances for the high outgassing cases than the model suggests (solid line). These simulations were made for the high ionisation rate expected close to the Sun, while the model was calculated for 1 AU conditions. A disagreement is therefore expected.

The minimum outgassing rate for which a bow shock forms (i.e. when $R_{BS}$ goes to zero) is about 3.5 × 10²⁷ s⁻¹ according to
<figure><img src="image_5.png" /><figcaption>Fig. 4. Bow shocks of comets for varying outgassing rates. For comparison, the bow shocks of Mars, Venus, Earth, Saturn, and Jupiter are indicated by dotted lines, and the approximate bow shock shapes of comets 1P, 19P, 21P and 26 are plotted with dashed lines.</figcaption></figure>

Eq. (3), shown in Fig. 2, for the conditions at 1 AU. For instance, decreasing the pick-up time $t_p$ by 20% lowers this to about $2 \times 10^{27}$ s$^{-1}$, and increasing by 20% increases this to about $5.5 \times 10^{27}$ s$^{-1}$. The infant bow shocks reported by Gunell et al. (2018) and Goetz et al. (2021) are mostly found below a few hundred kilometers.

Each simulation that is included in Fig. 3 has used a specific set of input parameters, such as solar wind density $n_{sw}$ and velocity $v_{sw}$, ionisation rates $v$, and the cometary neutral gas flow velocity $u_n$ for their models. Similarly, the spacecraft-measured crossings have all occurred during slightly different ambient condition. The standoff distance of all bow shocks might therefore potentially be normalised to common conditions representative of a heliocentric distance of 1 AU, but given the many parameters involved (see Eq. (3)), this would be rather non-trivial and is not done in this paper.

For a weakly outgassing comet, the atmosphere is optically thin, while for a more active comet, the atmosphere may become optically thick such that solar EUV will not reach and ionise the inner part of the coma (Bhardwaj 2003; Vigren & Galand 2013; Beth et al. 2019). However, this should still occur relatively close to the nucleus and should not impact the location of the bow shock, which is at significantly larger distances.

Furthermore, we provided a simple bow shock shape model, assuming that the bow shock takes the shape of a conic section, which seems to generally be the case for solar system bow shocks. The bow shock model is dependent on the outgassing rate, the ionisation rate, the flow velocity, and the upstream solar wind conditions. We assumed that the shape of the bow shock is rotationally symmetric, which might not be realistic as the bow shock location will vary with the upstream IMF and with the direction of the convective electric field, the solar wind flow shears, and the asymmetric neutral outgassing from an inhomogeneous surface of the comet nucleus, but also because of the asymmetry of the cometary ion pick-up process within the comet-induced magnetosphere. This still gives a good estimate of the bow shock shape, and it compares reasonably well to shapes presented in simulation papers.

Overall, our comparison between observations or simulations and the gas-dynamic model agree well. However, the uncertainty in the expected bow shock location is large and

depends on several factors, as mentioned above. Nevertheless, this model works over a range of magnitudes, which allows us to predict the bow shock location reasonably well. It may seem surprising that the same approximation holds for a broad range of gas-production rates and that no marked change occurs for the transitions from the fluid regime for a very active comet to the kinetic regime for a weakly outgassing comet. The reason is that the bow shock position is largely dictated by solar wind mass flux and cometary neutral outflow speed, but does not depend directly on cometary neutral density or on the corresponding ion density, which determines the plasma regime.

## 5. What to expect for cometary flybys

Figure 5 presents an example of the plasma density and the location of the bow shock that would be encountered during a cometary flyby for two cases of outgassing rates (a low rate of $10^{28}$ s$^{-1}$ and a high rate of $10^{30}$ s$^{-1}$). The plasma density was calculated using two different models. First, the plasma density was calculated from the ionisation of a Haser model of the neutral outgassing profile,

$$n_i(r) = \frac{vQ}{4\pi r u^2}, \qquad (7)$$

where $u \approx u_n$ is the radial ion flow velocity, as shown by Vigren et al. (e.g. 2015); Galand et al. (e.g. 2016). For higher outgassing rates, this model becomes less appropriate because it does not take electron-ion dissociative recombination into account, which significantly reduces the plasma density (Gombosi 2015; Beth et al. 2019). We therefore also included the model by Gombosi (2015),

$$n_{i,DR}(r) = \left( \sqrt{1 + \frac{Q}{Q_0}} - 1 \right) \frac{u_n}{2\alpha r}, \qquad (8)$$

where $Q_0 = \frac{\pi u_n^3}{v\alpha}$, and $\alpha$ is the water ion-electron recombination rate, for which we simply adopted the value $\alpha = 4.3 \sqrt{300}/T_e \times 10^{-13} \text{ m}^3 \text{ s}^{-1}$ (Beth et al. 2019) with $T_e = 150$ K. For the high outgassing case, the difference in density between the two models is almost an order of magnitude, while
<figure><img src="image_6.png" /><figcaption>Fig. 4. Bow shocks of comets for varying outgassing rates. For comparison, the bow shocks of Mars, Venus, Earth, Saturn, and Jupiter are indicated by dotted lines, and the approximate bow shock shapes of comets 1P, 19P, 21P and 26 are plotted with dashed lines.</figcaption></figure>

Eq. (3), shown in Fig. 2, for the conditions at 1 AU. For instance, decreasing the pick-up time $t_p$ by 20% lowers this to about $2 \times 10^{27}$ s$^{-1}$, and increasing by 20% increases this to about $5.5 \times 10^{27}$ s$^{-1}$. The infant bow shocks reported by Gunell et al. (2018) and Goetz et al. (2021) are mostly found below a few hundred kilometers.

Each simulation that is included in Fig. 3 has used a specific set of input parameters, such as solar wind density $n_{sw}$ and velocity $v_{sw}$, ionisation rates $v$, and the cometary neutral gas flow velocity $u_n$ for their models. Similarly, the spacecraft-measured crossings have all occurred during slightly different ambient condition. The standoff distance of all bow shocks might therefore potentially be normalised to common conditions representative of a heliocentric distance of 1 AU, but given the many parameters involved (see Eq. (3)), this would be rather non-trivial and is not done in this paper.

For a weakly outgassing comet, the atmosphere is optically thin, while for a more active comet, the atmosphere may become optically thick such that solar EUV will not reach and ionise the inner part of the coma (Bhardwaj 2003; Vigren & Galand 2013; Beth et al. 2019). However, this should still occur relatively close to the nucleus and should not impact the location of the bow shock, which is at significantly larger distances.

Furthermore, we provided a simple bow shock shape model, assuming that the bow shock takes the shape of a conic section, which seems to generally be the case for solar system bow shocks. The bow shock model is dependent on the outgassing rate, the ionisation rate, the flow velocity, and the upstream solar wind conditions. We assumed that the shape of the bow shock is rotationally symmetric, which might not be realistic as the bow shock location will vary with the upstream IMF and with the direction of the convective electric field, the solar wind flow shears, and the asymmetric neutral outgassing from an inhomogeneous surface of the comet nucleus, but also because of the asymmetry of the cometary ion pick-up process within the comet-induced magnetosphere. This still gives a good estimate of the bow shock shape, and it compares reasonably well to shapes presented in simulation papers.

Overall, our comparison between observations or simulations and the gas-dynamic model agree well. However, the uncertainty in the expected bow shock location is large and

depends on several factors, as mentioned above. Nevertheless, this model works over a range of magnitudes, which allows us to predict the bow shock location reasonably well. It may seem surprising that the same approximation holds for a broad range of gas-production rates and that no marked change occurs for the transitions from the fluid regime for a very active comet to the kinetic regime for a weakly outgassing comet. The reason is that the bow shock position is largely dictated by solar wind mass flux and cometary neutral outflow speed, but does not depend directly on cometary neutral density or on the corresponding ion density, which determines the plasma regime.

## 5. What to expect for cometary flybys

Figure 5 presents an example of the plasma density and the location of the bow shock that would be encountered during a cometary flyby for two cases of outgassing rates (a low rate of $10^{28}$ s$^{-1}$ and a high rate of $10^{30}$ s$^{-1}$). The plasma density was calculated using two different models. First, the plasma density was calculated from the ionisation of a Haser model of the neutral outgassing profile,

$$n_i(r) = \frac{vQ}{4\pi r u^2}, \qquad (7)$$

where $u \approx u_n$ is the radial ion flow velocity, as shown by Vigren et al. (e.g. 2015); Galand et al. (e.g. 2016). For higher outgassing rates, this model becomes less appropriate because it does not take electron-ion dissociative recombination into account, which significantly reduces the plasma density (Gombosi 2015; Beth et al. 2019). We therefore also included the model by Gombosi (2015),

$$n_{i,DR}(r) = \left( \sqrt{1 + \frac{Q}{Q_0}} - 1 \right) \frac{u_n}{2\alpha r}, \qquad (8)$$

where $Q_0 = \frac{\pi u_n^3}{v\alpha}$, and $\alpha$ is the water ion-electron recombination rate, for which we simply adopted the value $\alpha = 4.3 \sqrt{300}/T_e \times 10^{-13} \text{ m}^3 \text{ s}^{-1}$ (Beth et al. 2019) with $T_e = 150$ K. For the high outgassing case, the difference in density between the two models is almost an order of magnitude, while
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Ionospheric density as a function of cometocentric distance from two models (solid and dashed lines, see text) and assuming that the density starts to decrease as $1/r^2$ at $10^4$ km. The density is calculated for two outgassing rates (red and blue). The bow shock distance at the terminator (which is roughly twice the subsolar standoff distance) for the two outgassing rates is indicated (vertical dash-dotted lines). The grey shaded region indicates the solar wind density at 1 AU. The ion density measured by Giotto at comet 1P/Halley is also included (green dots).</figcaption></figure>

for the low outgassing rate, the difference is less significant. The measured ion density from 1P/Halley (outgassing rate of ~$10^{30}$ s⁻¹; Altwegg et al. 1993) is also included in Fig. 5. This agrees well with $n_{i,DR}$, at least within a cometocentric distance of $10^4$ km. Inside of this point, the density falls off as $1/r$, which has also been observed with Rosetta at comet 67P (Edberg et al. 2015, 2022). Beyond this, the plasma density increases locally before it again decreases with a steeper $1/r^2$ trend.

Beyond ~$10^4$ km, the ions and neutrals should be dynamically decoupled, and the plasma velocity at 1P/Halley was indeed measured to increase linearly with $r$ beyond this range, out to about $3 \times 10^4$ km (Schwenn et al. 1987; Altwegg et al. 1993). This resulted in a more rapid decrease in the density with distance. The $1/r^2$ trend is indicated by dotted lines and follows the Giotto data quite well for the high outgassing case. For the low outgassing rate, we assumed that this steeper trend also starts at $10^4$ km, although the location of the trend break likely depends on the outgassing rate, is therefore possibly closer to the nucleus for lower outgassing, and needs not necessarily be abrupt.

As the cometocentric distance increases, the cometary plasma density eventually decreases to the point at which it becomes equal to the solar wind density, about $5 \text{ cm}^{-3}$ (illustrated as the grey shaded region in Fig. 5). The cometocentric distance of this point can be compared to the expected bow shock location for the two outgassing rates, as calculated using Eq. (3). The bow shock distance at the terminator (when $\theta = 90^\circ$) is shown by the two vertical dash-dotted lines. For the low outgassing case, the bow shock is well within the distance within which the cometary plasma density has dropped to solar wind values, but the $1/r^2$ trend might start closer to the nucleus than $10^4$ km. For the high outgassing case, the bow shock is instead well outside of the point where $n_i = n_{sw}$ when we assume that the $1/r^2$ trend is present.

For a spacecraft flyby speed of 10 km s⁻¹ or 70 km s⁻¹ past a comet with outgassing rate of $10^{28}$ s⁻¹ or $10^{30}$ s⁻¹, the time spent within the bow shock of a comet can vary quite significantly, as can be read off from the lower two horizontal axes

in Fig. 5. For a fast flyby past a weakly outgassing comet in the terminator plane and if there is no significant aberration, the time moving from the bow shock to closest approach and out to the bow shock on the far side is about 10 minutes. For the combination of a slow flyby and a very active comet, the spacecraft could spend up to 3 days within the bow shock. The most likely flyby velocity for Comet Interceptor is predicted (given the statistics of currently known comets) to be 50 km s⁻¹, and combined with an outgassing rate of $10^{28}$ s⁻¹, the time spent within the bow shock would be approximately 13 min.

Beyond the bow shock, there can still be plasma signatures of the comet in terms of pick-up ions and associated plasma waves. From measurements at 1P/Halley, waves like this were observed at approximately twice the bow shock stand-off distance (Glassmeier et al. 1989); at comet 26P/Grigg-Skjellerup, solar wind contaminated by pick-up ions was found already at a distance of $6 \times 10^5$ km, that is, more than an order of magnitude farther out than the bow shock (Glassmeier & Neubauer 1993), and at 19P/Borrelly, pick-up ions were also observed ~$6 \times 10^5$ km upstream, that is, about an order of magnitude farther out than the bow shock (Young et al. 2004; Delamere 2006).

## 6. Conclusions

Our compiled data set of spacecraft measurements and simulation results of cometary bow shock locations has allowed us to study the dependence of the bow shock size on the comet outgassing rate. The good agreement in the bow shock stand-off distance between the data and the model developed by Koenders et al. (2013), based on Galeev et al. (1985) and Biermann et al. (1967; Eq. (3)), provides confidence that the model is sufficiently good for these purposes.

The scale size of comets is found to vary by four orders of magnitude, from ~$10^2$ km to $10^6$ km, for outgassing rates in the range of $10^{27}-10^{31}$ s⁻¹ and an ionisation rate, flow velocity, and upstream solar wind conditions typical of those at 1 AU.
<figure><img src="image_7.png" /><figcaption>**Fig. 5.** Ionospheric density as a function of cometocentric distance from two models (solid and dashed lines, see text) and assuming that the density starts to decrease as $1/r^2$ at $10^4$ km. The density is calculated for two outgassing rates (red and blue). The bow shock distance at the terminator (which is roughly twice the subsolar standoff distance) for the two outgassing rates is indicated (vertical dash-dotted lines). The grey shaded region indicates the solar wind density at 1 AU. The ion density measured by Giotto at comet 1P/Halley is also included (green dots).</figcaption></figure>

for the low outgassing rate, the difference is less significant. The measured ion density from 1P/Halley (outgassing rate of ~$10^{30}$ s⁻¹; Altwegg et al. 1993) is also included in Fig. 5. This agrees well with $n_{i,DR}$, at least within a cometocentric distance of $10^4$ km. Inside of this point, the density falls off as $1/r$, which has also been observed with Rosetta at comet 67P (Edberg et al. 2015, 2022). Beyond this, the plasma density increases locally before it again decreases with a steeper $1/r^2$ trend.

Beyond ~$10^4$ km, the ions and neutrals should be dynamically decoupled, and the plasma velocity at 1P/Halley was indeed measured to increase linearly with $r$ beyond this range, out to about $3 \times 10^4$ km (Schwenn et al. 1987; Altwegg et al. 1993). This resulted in a more rapid decrease in the density with distance. The $1/r^2$ trend is indicated by dotted lines and follows the Giotto data quite well for the high outgassing case. For the low outgassing rate, we assumed that this steeper trend also starts at $10^4$ km, although the location of the trend break likely depends on the outgassing rate, is therefore possibly closer to the nucleus for lower outgassing, and needs not necessarily be abrupt.

As the cometocentric distance increases, the cometary plasma density eventually decreases to the point at which it becomes equal to the solar wind density, about $5 \text{ cm}^{-3}$ (illustrated as the grey shaded region in Fig. 5). The cometocentric distance of this point can be compared to the expected bow shock location for the two outgassing rates, as calculated using Eq. (3). The bow shock distance at the terminator (when $\theta = 90^\circ$) is shown by the two vertical dash-dotted lines. For the low outgassing case, the bow shock is well within the distance within which the cometary plasma density has dropped to solar wind values, but the $1/r^2$ trend might start closer to the nucleus than $10^4$ km. For the high outgassing case, the bow shock is instead well outside of the point where $n_i = n_{sw}$ when we assume that the $1/r^2$ trend is present.

For a spacecraft flyby speed of 10 km s⁻¹ or 70 km s⁻¹ past a comet with outgassing rate of $10^{28}$ s⁻¹ or $10^{30}$ s⁻¹, the time spent within the bow shock of a comet can vary quite significantly, as can be read off from the lower two horizontal axes

in Fig. 5. For a fast flyby past a weakly outgassing comet in the terminator plane and if there is no significant aberration, the time moving from the bow shock to closest approach and out to the bow shock on the far side is about 10 minutes. For the combination of a slow flyby and a very active comet, the spacecraft could spend up to 3 days within the bow shock. The most likely flyby velocity for Comet Interceptor is predicted (given the statistics of currently known comets) to be 50 km s⁻¹, and combined with an outgassing rate of $10^{28}$ s⁻¹, the time spent within the bow shock would be approximately 13 min.

Beyond the bow shock, there can still be plasma signatures of the comet in terms of pick-up ions and associated plasma waves. From measurements at 1P/Halley, waves like this were observed at approximately twice the bow shock stand-off distance (Glassmeier et al. 1989); at comet 26P/Grigg-Skjellerup, solar wind contaminated by pick-up ions was found already at a distance of $6 \times 10^5$ km, that is, more than an order of magnitude farther out than the bow shock (Glassmeier & Neubauer 1993), and at 19P/Borrelly, pick-up ions were also observed ~$6 \times 10^5$ km upstream, that is, about an order of magnitude farther out than the bow shock (Young et al. 2004; Delamere 2006).

## 6. Conclusions

Our compiled data set of spacecraft measurements and simulation results of cometary bow shock locations has allowed us to study the dependence of the bow shock size on the comet outgassing rate. The good agreement in the bow shock stand-off distance between the data and the model developed by Koenders et al. (2013), based on Galeev et al. (1985) and Biermann et al. (1967; Eq. (3)), provides confidence that the model is sufficiently good for these purposes.

The scale size of comets is found to vary by four orders of magnitude, from ~$10^2$ km to $10^6$ km, for outgassing rates in the range of $10^{27}-10^{31}$ s⁻¹ and an ionisation rate, flow velocity, and upstream solar wind conditions typical of those at 1 AU.
Although the bow shock location depends on many different parameters, the outgassing rate varies most strongly for the different comets.

A spacecraft flyby through a cometary plasma environment, such as the planned Comet Interceptor mission, could last from anywhere between tens of minutes to days within the bow shock, depending on the comet target outgassing rate, ionisation rate, flow speed of the neutrals and ions, and the relative flyby speed. To ensure as much plasma measurements as possible, the selected target comet should therefore preferably be relatively active, and the flyby velocity should be as low as possible.

For the practical purpose of planning the flyby timeline, there can still be phenomena (pick-up ions and plasma waves) upstream of the bow shock that are interesting to study and therefore should not be ignored when planning the operations. Because of the uncertainties on the model, the solar wind conditions, and the comet activity, it is of uttermost importance to allow for a sizeable margin. Measurements should start well ahead of the expected bow shock crossing to ensure that both upstream phenomena and the transition into the comet-dominated plasma regime are captured.

*Acknowledgements*. The work of HN was supported by SNSA grant 132/19. The work of NE was supported by SNSA grant 2021-00047 and VR grant 2020-03962. JDK acknowledges support from BELSPO via ESA/PRODEX PEA 4000139830.

## References

Alho, M., Jarvinen, R., Wedlund, C. S., et al. 2021, MNRAS, 506, 4735

Altwegg, K., Balsiger, H., Geiss, J., et al. 1993, A&A, 279, 260

Bame, S. J., Anderson, R. C., Asbridge, J. R., et al. 1986, Science, 232, 356

Behar, E., & Henri, P. 2023, A&A, 671, A144

Behar, E., Fatemi, S., Henri, P., & Holmström, M. 2022, Ann. Geophys., 40, 281

Benna, M., & Mahaffy, P. R. 2006, Geophys. Res. Lett., 33

Benna, M., & Mahaffy, P. 2007, Planet. Space Sci., 55, 1031

Beth, A., Galand, M., & Heritier, K. L. 2019, A&A, 630, A47

Bhardwaj, A. 2003, Geophys. Res. Lett., 30, 24

Biermann, L., Brosowski, B., & Schmidt, H. 1967, Solar Phys., 1, 254

Coates, A. J., Mazelle, C., & Neubauer, F. M. 1997, J. Geophys. Res. Space Phys., 102, 7105

Deca, J., Divin, A., Henri, P., et al. 2017, Phys. Rev. Lett., 118, 205101

Deca, J., Henri, P., Divin, A., et al. 2019, Phys. Rev. Lett., 123, 055101

Delamere, P. A. 2006, J. Geophys. Res., 111, 12217

Edberg, N. J. T., Lester, M., Cowley, S. W. H., & Eriksson, A. I. 2008, J. Geophys. Res., 113, A08206

Edberg, N. J. T., Andrews, D. J., Bertucci, C., et al. 2015, J. Geophys. Res. Space Phys., 120, 8884, 2015JA021373

Edberg, N. J. T., Johansson, F. L., Eriksson, A. I., et al. 2022, A&A, 663, A42

Farris, M. H., Petrinec, S. M., & Russell, C. T. 1991, Geophys. Res. Lett., 18, 1821

Flammer, K. R., & Mendis, D. A. 1993, J. Geophys. Res., 98, 21003

Galand, M., Héritier, K. L., Odelstad, E., et al. 2016, MNRAS, 462, S331

Galeev, A. A. 1985, Adv. Space Res., 5, 155

Galeev, A. A., Cravens, T. E., & Gombosi, T. I. 1985, ApJ, 289, 807

Glassmeier, K.-H., & Neubauer, F. M. 1993, J. Geophys. Res., 98, 20921

Glassmeier, K.-H., Coates, A. J., Aculqa, M. H., et al. 1989, J. Geophys. Res., 94, 37

Glassmeier, K.-H., Boehnhardt, H., Koschny, D., Kührt, E., & Richter, I. 2007, Space Sci. Rev., 128, 1

Goetz, C., Gunell, H., Johansson, F., et al. 2021, Ann. Geophys., 39, 379

Goetz, C., Behar, E., Beth, A., et al. 2022, Space Sci. Rev., 218, 65

Gombosi, T. I. 2015, Geophysical Monograph Series, 207, 169

Gringauz, K. I., Gombosi, T. I., Remizov, A. P., et al. 1986, Nature, 321, 282

Gunell, H., Goetz, Charlotte, Wedlund, Cyril Simon, et al. 2018, A&A, 619, L2

Hansen, K., Bagdonat, T., Motschmann, U., et al. 2007, Space Sci. Rev., 128, 133

Hirao, K., & Itoh, T. 1986, Nature, 321, 294

Huang, Z., Tóth, G., Gombosi, T. I., et al. 2016, J. Geophys. Res. Space Phys., 121, 4247

Huddleston, D. E., Russell, C. T., Kivelson, M. G., Khurana, K. K., & Bennett, L. 1998, J. Geophys. Res. Planets, 103, 20075

Jia, Y. D., Combi, M. R., Hansen, K. C., et al. 2008, Icarus, 196, 249

Jia, Y. D., Russell, C. T., Jian, L. K., et al. 2009, ApJ, 696

Jia, Y. D., Russell, C. T., Liu, W., & Shou, Y. S. 2014, ApJ, 796

Johnstone, A. D., Coates, A. J., Huddleston, D. E., et al. 1993, A&A, 273, L1

Jones, D. E., Smith, E. J., Slavin, J. A., et al. 1986, Geophys. Res. Lett., 13, 243

Jones, G., Snodgrass, C., & Tubiana, C. 2022, Space Sci. Rev., submitted

Koenders, C., Glassmeier, K.-H., Richter, I., Motschmann, U., & Rubin, M. 2013, Planet. Space Sci., 87, 85

Koenders, C., Glassmeier, K.-H., Richter, I., Ranocha, H., & Motschmann, U. 2015, Planet. Space Sci., 105, 101

Lindkvist, J., Hamrin, M., Gunell, H., et al. 2018, A&A, 616, A81

Martinecz, C., Boesswetter, A., Fränz, M., et al. 2009, J. Geophys. Res. Planets, 114

Masters, A., Achilleos, N., Dougherty, M. K., et al. 2008, J. Geophys. Res., 113, A10210

Mazelle, C., Winterhalter, D., Sauer, K., et al. 2004, Space Sci. Rev., 111, 115

Mukai, T., Miyake, W., Terasawa, T., Kitayama, M., & Hirao, K. 1986, Geophys. Res. Lett., 13, 829

Reinhard, R. 1986, Nature, 321, 313

Richter, I., Koenders, C., Glassmeier, K. H., Tsurutani, B. T., & Goldstein, R. 2011, Planet. Space Sci., 59, 691

Rosenvinge, T. T. V., Brandt, J. C., & Farquhar, R. W. 1986, Science, 232, 353

Rubin, M., Hansen, K. C., Gombosi, T. I., et al. 2009, Icarus, 199, 505

Sagdeev, R. Z., Blamont, J. E., Galeev, A. A., et al. 1986, Nature, 321, 259

Schleicher, D. G., Woodney, L. M., & Millis, R. L. 2003, Icarus, 162, 415

Schwenn, R., Ip, W. H., Rosenbauer, H., et al. 1987, A&A, 187, 160

Shou, Y., Combi, M., Jia, Y.-D., et al. 2015, ApJ, 809, 156

Simon Wedlund, C., Alho, M., Gronoff, G., et al. 2017, A&A, 604, A73

Snodgrass, C., & Jones, G. H. 2019, Nat. Commun., 10, 5418

Vigren, E., & Galand, M. 2013, ApJ, 772, 33

Vigren, E., Galand, M., Eriksson, A. I., et al. 2015, ApJ, 812, 9

Young, D. T., Crary, F. J., Nordholt, J. E., et al. 2004, Icarus, 167, 80
Although the bow shock location depends on many different parameters, the outgassing rate varies most strongly for the different comets.

A spacecraft flyby through a cometary plasma environment, such as the planned Comet Interceptor mission, could last from anywhere between tens of minutes to days within the bow shock, depending on the comet target outgassing rate, ionisation rate, flow speed of the neutrals and ions, and the relative flyby speed. To ensure as much plasma measurements as possible, the selected target comet should therefore preferably be relatively active, and the flyby velocity should be as low as possible.

For the practical purpose of planning the flyby timeline, there can still be phenomena (pick-up ions and plasma waves) upstream of the bow shock that are interesting to study and therefore should not be ignored when planning the operations. Because of the uncertainties on the model, the solar wind conditions, and the comet activity, it is of uttermost importance to allow for a sizeable margin. Measurements should start well ahead of the expected bow shock crossing to ensure that both upstream phenomena and the transition into the comet-dominated plasma regime are captured.

*Acknowledgements*. The work of HN was supported by SNSA grant 132/19. The work of NE was supported by SNSA grant 2021-00047 and VR grant 2020-03962. JDK acknowledges support from BELSPO via ESA/PRODEX PEA 4000139830.

## References

Alho, M., Jarvinen, R., Wedlund, C. S., et al. 2021, MNRAS, 506, 4735

Altwegg, K., Balsiger, H., Geiss, J., et al. 1993, A&A, 279, 260

Bame, S. J., Anderson, R. C., Asbridge, J. R., et al. 1986, Science, 232, 356

Behar, E., & Henri, P. 2023, A&A, 671, A144

Behar, E., Fatemi, S., Henri, P., & Holmström, M. 2022, Ann. Geophys., 40, 281

Benna, M., & Mahaffy, P. R. 2006, Geophys. Res. Lett., 33

Benna, M., & Mahaffy, P. 2007, Planet. Space Sci., 55, 1031

Beth, A., Galand, M., & Heritier, K. L. 2019, A&A, 630, A47

Bhardwaj, A. 2003, Geophys. Res. Lett., 30, 24

Biermann, L., Brosowski, B., & Schmidt, H. 1967, Solar Phys., 1, 254

Coates, A. J., Mazelle, C., & Neubauer, F. M. 1997, J. Geophys. Res. Space Phys., 102, 7105

Deca, J., Divin, A., Henri, P., et al. 2017, Phys. Rev. Lett., 118, 205101

Deca, J., Henri, P., Divin, A., et al. 2019, Phys. Rev. Lett., 123, 055101

Delamere, P. A. 2006, J. Geophys. Res., 111, 12217

Edberg, N. J. T., Lester, M., Cowley, S. W. H., & Eriksson, A. I. 2008, J. Geophys. Res., 113, A08206

Edberg, N. J. T., Andrews, D. J., Bertucci, C., et al. 2015, J. Geophys. Res. Space Phys., 120, 8884, 2015JA021373

Edberg, N. J. T., Johansson, F. L., Eriksson, A. I., et al. 2022, A&A, 663, A42

Farris, M. H., Petrinec, S. M., & Russell, C. T. 1991, Geophys. Res. Lett., 18, 1821

Flammer, K. R., & Mendis, D. A. 1993, J. Geophys. Res., 98, 21003

Galand, M., Héritier, K. L., Odelstad, E., et al. 2016, MNRAS, 462, S331

Galeev, A. A. 1985, Adv. Space Res., 5, 155

Galeev, A. A., Cravens, T. E., & Gombosi, T. I. 1985, ApJ, 289, 807

Glassmeier, K.-H., & Neubauer, F. M. 1993, J. Geophys. Res., 98, 20921

Glassmeier, K.-H., Coates, A. J., Aculqa, M. H., et al. 1989, J. Geophys. Res., 94, 37

Glassmeier, K.-H., Boehnhardt, H., Koschny, D., Kührt, E., & Richter, I. 2007, Space Sci. Rev., 128, 1

Goetz, C., Gunell, H., Johansson, F., et al. 2021, Ann. Geophys., 39, 379

Goetz, C., Behar, E., Beth, A., et al. 2022, Space Sci. Rev., 218, 65

Gombosi, T. I. 2015, Geophysical Monograph Series, 207, 169

Gringauz, K. I., Gombosi, T. I., Remizov, A. P., et al. 1986, Nature, 321, 282

Gunell, H., Goetz, Charlotte, Wedlund, Cyril Simon, et al. 2018, A&A, 619, L2

Hansen, K., Bagdonat, T., Motschmann, U., et al. 2007, Space Sci. Rev., 128, 133

Hirao, K., & Itoh, T. 1986, Nature, 321, 294

Huang, Z., Tóth, G., Gombosi, T. I., et al. 2016, J. Geophys. Res. Space Phys., 121, 4247

Huddleston, D. E., Russell, C. T., Kivelson, M. G., Khurana, K. K., & Bennett, L. 1998, J. Geophys. Res. Planets, 103, 20075

Jia, Y. D., Combi, M. R., Hansen, K. C., et al. 2008, Icarus, 196, 249

Jia, Y. D., Russell, C. T., Jian, L. K., et al. 2009, ApJ, 696

Jia, Y. D., Russell, C. T., Liu, W., & Shou, Y. S. 2014, ApJ, 796

Johnstone, A. D., Coates, A. J., Huddleston, D. E., et al. 1993, A&A, 273, L1

Jones, D. E., Smith, E. J., Slavin, J. A., et al. 1986, Geophys. Res. Lett., 13, 243

Jones, G., Snodgrass, C., & Tubiana, C. 2022, Space Sci. Rev., submitted

Koenders, C., Glassmeier, K.-H., Richter, I., Motschmann, U., & Rubin, M. 2013, Planet. Space Sci., 87, 85

Koenders, C., Glassmeier, K.-H., Richter, I., Ranocha, H., & Motschmann, U. 2015, Planet. Space Sci., 105, 101

Lindkvist, J., Hamrin, M., Gunell, H., et al. 2018, A&A, 616, A81

Martinecz, C., Boesswetter, A., Fränz, M., et al. 2009, J. Geophys. Res. Planets, 114

Masters, A., Achilleos, N., Dougherty, M. K., et al. 2008, J. Geophys. Res., 113, A10210

Mazelle, C., Winterhalter, D., Sauer, K., et al. 2004, Space Sci. Rev., 111, 115

Mukai, T., Miyake, W., Terasawa, T., Kitayama, M., & Hirao, K. 1986, Geophys. Res. Lett., 13, 829

Reinhard, R. 1986, Nature, 321, 313

Richter, I., Koenders, C., Glassmeier, K. H., Tsurutani, B. T., & Goldstein, R. 2011, Planet. Space Sci., 59, 691

Rosenvinge, T. T. V., Brandt, J. C., & Farquhar, R. W. 1986, Science, 232, 353

Rubin, M., Hansen, K. C., Gombosi, T. I., et al. 2009, Icarus, 199, 505

Sagdeev, R. Z., Blamont, J. E., Galeev, A. A., et al. 1986, Nature, 321, 259

Schleicher, D. G., Woodney, L. M., & Millis, R. L. 2003, Icarus, 162, 415

Schwenn, R., Ip, W. H., Rosenbauer, H., et al. 1987, A&A, 187, 160

Shou, Y., Combi, M., Jia, Y.-D., et al. 2015, ApJ, 809, 156

Simon Wedlund, C., Alho, M., Gronoff, G., et al. 2017, A&A, 604, A73

Snodgrass, C., & Jones, G. H. 2019, Nat. Commun., 10, 5418

Vigren, E., & Galand, M. 2013, ApJ, 772, 33

Vigren, E., Galand, M., Eriksson, A. I., et al. 2015, ApJ, 812, 9

Young, D. T., Crary, F. J., Nordholt, J. E., et al. 2004, Icarus, 167, 80