# A Low-Cost, Easy-to-Manufacture, Flexible, Multi-Taxel Tactile Sensor
and its Application to In-Hand Object Recognition

Tessa J. Pannen* Steffen Puhlmann* Oliver Brock  *These authors contributed
equally to this work.
All authors are with the Robotics and Biology Laboratory, Technische
Universität Berlin, Germany. We gratefully acknowledge financial support by
the European Commission (SOMA, H2020-ICT-645599), the Deutsche
Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s
Excellence Strategy - EXC 2002/1 ”Science of Intelligence” - project number
390523135 and German Priority Program DFG-SPP 2100 “Soft Material Robotic
Systems” - project number 405033880.

###### Abstract

Soft robotics is an emerging field that yields promising results for tasks
that require safe and robust interactions with the environment or with humans,
such as grasping, manipulation, and human-robot interaction. Soft robots rely
on intrinsically compliant components and are difficult to equip with
traditional, rigid sensors which would interfere with their compliance. We
propose a highly flexible tactile sensor that is low-cost and easy to
manufacture while measuring contact pressures independently from 14 taxels.
The sensor is built from piezoresistive fabric for highly sensitive,
continuous responses and from a custom-designed flexible printed circuit board
which provides a high taxel density. From these taxels, location and intensity
of contact with the sensor can be inferred. In this paper, we explain the
design and manufacturing of the proposed sensor, characterize its input-output
relation, evaluate its effects on compliance when equipped to the silicone-
based pneumatic actuators of the soft robotic RBO Hand 2, and demonstrate that
the sensor provides rich and useful feedback for learning-based in-hand object
recognition.

##  I Introduction

Soft robotic hands demonstrate robust grasping and manipulation capabilities
thanks to their intrinsic compliance [1, 2, 3]. However, sensorization of
these soft hands is challenging, as suitable sensors must provide sensing
abilities without detrimental effects on compliance. Soft sensors therefore
need to be based on flexible materials, however, useful designs that provide
rich feedback often rely on complex fabrication or costly components.

We propose a novel, highly compliant tactile sensor (Fig. 1) suitable for soft
robotics that is low-cost and very easy to manufacture while providing rich
sensory feedback with high spatial resolution, enabling learning-based object
recognition. The sensor relies on piezoresistive fabric which changes its
electrical resistance under contact pressure. It is highly sensitive while
measuring a wide range of contact pressures. Its sensitivity and measuring
range can be adjusted to support a large variety of applications with
different contact intensities. To effectively distinguish contact patterns,
our sensor possesses 14 taxels realized by a custom-designed flexible printed
circuit board (flexPCB). This flexPCB tightly integrates electrodes and wiring
in a thin layout.

We evaluate the sensor in grasping experiments with the highly compliant,
anthropomorphic RBO Hand 2 [4] and show that it provides useful tactile
information, enabling classification of grasped objects. We also show that the
sensor has only limited effect on the hand’s compliance, so that passive shape
adaptation during grasping is still attainable. Finally, we compare object
classification based on the tactile sensor to classification based on air
pressure sensing in the hand’s pneumatic actuators. Our results show that
thanks to its superior spatial resolution, our sensor outperforms pressure-
based object classification, highlighting that our cost-effective and simple
design yields rich feedback.

Figure 1: A flexible tactile sensor: palmar side shows the grounding electrode
made of Medtex P130 fabric (top), dorsal side shows the flexPCB integrating
electrodes and wiring (bottom). The sensor has 14 independent taxels (the most
proximal electrode/taxel is not used).

##  II Related Work

Related works presented a wide range of tactile sensors based on a variety of
technologies [5, 6, 7]. We now discuss promising research that contributes to
achieving dense tactile sensing while limiting detrimental effects on the
compliance of a sensorized soft actuator.

Microchannels filled with liquid metal change their resistance under
deformation, enabling the measurement of multi-axis strain and pressure [8,
9]. Several such sensors can be combined to infer contact information on the
entire actuator, including deformation, forces, and contact location [10, 11].
These sensors have only minor effect on the compliance of a sensorized
actuator but they have not yet been shown to provide dense tactile sensing.

Acoustic sensing relies on sound propagating through a soft continuum actuator
to infer information about contact location, contact intensity, and contacted
material [12, 13]. Since the required microphone and optional speaker can be
embedded into the cavity of the pneumatic actuator, the effect on compliance
is negligible. However, only single contact locations can be measured. The
sensor presented in this paper allows concurrent sensing of many contacts.

Optical tactile sensors use cameras to record deformation-based changes of
marker patterns [14, 15] or of the surface relief [16] on the inside of soft
actuators. Embedding these cameras requires a rigid base for the soft
actuator, limiting the compliance and design choices.

Elastomers enriched with carbon black particles change their resistance under
deformation [17, 18]. However, these particles can affect the elastomer’s
mechanical properties, so that increasing sensitivity comes at the expense of
lower flexibility [18].

A different design for soft tactile sensors sandwiches piezoresistive foil or
fabric between electrode layers. Individual patches of the piezoresistive
material encode the taxel structure [19, 20]. Placed on a soft robotic finger,
these sensors can infer gasping forces and locations. However, they have
relatively low tactile density and published designs often have external cable
routing, possibly complicating operation in proximity to obstacles.

Alternatively, the electrode layers encode the taxel structure [21, 22]. The
fabrication of these sensors is often complex, for example requiring etching
of the electrode layers. However, these designs offer relatively high taxel
density while promising not to reduce compliance substantially. We therefore
base our work on this kind of design. Our sensor, described in the next
section, offers simple manufacturing and dense tactile sensing, overcoming the
problems of existing tactile sensors based on this design.

##  III Tactile Sensor Design and Manufacturing

Figure 2: Fabrication of the tactile sensor: a custom-designed flexPCB
(bottom) which integrates fourteen independent taxels and wiring, a laser-cut
sheet of EeonTex piezoresistive fabric (center) which changes its electrical
resistance under contact pressure, and a laser-cut sheet of silver-coated
nylon fabric (top) which serves as a grounding electrode are connected by
laser-cut double-sided adhesive tape.

Our tactile sensor should have only little effect on the robot’s morphology
and compliance. It should exhibit a high spatial resolution and provide
continuous responses for a large range of contact pressures with high
sensitivity. At the same time, it should be based on low-cost materials and
very simple to manufacture. We now describe how we achieve these design goals.

The sensor is based on a piezoresistive material which is sandwiched between
two electrode layers for reading-out its electrical resistance (Fig. 1). We
divide the sensor into different regions (taxels) for which we determine the
resistance independently (electrodes of taxels are connected in parallel). To
minimize the sensor’s effect on morphology and compliance, we choose thin and
flexible materials, including conductive fabrics and a flexPCB.

The sensor design is based on EeonTex, a commercially available piezoresistive
fabric which is only $0.4\text{\,}\mathrm{m}\mathrm{m}$ thick and made of
nylon fabric coated with conductive polymers. EeonTex enables highly sensitive
responses, changing its electrical resistance by more than 99% under contact
pressure. This resistance is measured at different regions via multiple
electrodes (taxels) on a flexPCB and a single grounding electrode based on
Medtex P130, conductive fabric made of silver-plated nylon (Fig. 2).

The flexPCB provides a high spatial resolution of up to
$4.5\text{\,}\mathrm{m}\mathrm{m}$ at the fingertip, realized by 14 taxels
while integrating electrodes and cabling in a custom-designed, easily
accessible and low-cost component. At the fingertip, it has nine square
electrodes ($3\text{\,}\mathrm{m}\mathrm{m}$ side length) and five rectangular
taxels (3 $\times$ $13\text{\,}\mathrm{m}\mathrm{m}$) located proximally (Fig.
2). The flexPCB enables a clear and clutter-free design by keeping the sensing
surface void of solder joints and shifting cabling connections out of the
actuator’s workspace. These design choices result in a thin, compliant, and
highly sensitive sensor whose thickness ranges between
$1.0\text{\,}\mathrm{m}\mathrm{m}$ and $1.2\text{\,}\mathrm{m}\mathrm{m}$ due
to small variations in manufacturing.

Manufacturing the sensor is very simple and fast, allowing for rapid
prototyping and easy replacements. The fabrics and the flexPCB are connected
with double-sided adhesive tape (Fig. 2). The fabrics and the tape are laser-
cut with cut-outs in the tape that correspond to the shapes and locations of
the electrodes on the flexPCB. Finally, we solder cables to wire outlets at
the bottom of the flexPCB. Following this procedure, a sensor can be assembled
within minutes. The cost of the sensor is dominated by the commercially
manufactured flexPCB while the cost of a single piece drops substantially when
ordering large quantities. In our case, a batch size of 50 pieces resulted in
a sensor costing less than US$ 5, demonstrating the cost-effectiveness of our
design.

The sensor is glued to the palmar side of the silicone-based finger of the RBO
Hand 2 using Sil-Poxy silicone adhesive, with the flexPCB facing the actuator.
Since this wet adhesive would soak into the fabrics, we cover the sensor with
a thin, water resistant sleeve made of Dragon Skin 10 silicone.

##  IV Sensor Characterization

Figure 3: The sensor response can be adjusted by the choice of pull-up
resistance, trading measuring range with sensitivity (left). Sensor shows
repeatable hysteresis due to relaxation time of the piezoresistive fabric
(right). Solid lines show mean responses over ten trials, intervals indicate
standard deviation.

We now characterize the sensor’s response behavior and explain how it can be
adjusted to match task requirements, before we analyze the sensor’s effect on
compliance.

###  IV-A Response Behavior

We probe the sensor’s response behavior to analyze its sensitivity, measuring
range, and hysteresis. For this, we apply normal forces to a fingertip taxel,
using a square probe-tip of $3\text{\,}\mathrm{m}\mathrm{m}$ side length which
is attached to a force-torque sensor. Contact forces and sensor responses are
recorded at $20\text{\,}\mathrm{H}\mathrm{z}$, with each recording lasting
$3.5\text{\,}\mathrm{s}$.

In response to contact, the electrical resistance of the EeonTex fabric
changes at the location of the taxel. Instead of measuring this resistance
directly, we use a voltage divider and infer the response indirectly by
measuring the voltage between the conductive Medtex P30 fabric and the
electrode on the flexPCB that corresponds to this taxel. The relationship
between contact pressure and output voltage is described by the characteristic
curve. The sensor is most sensitive at the center of its measuring range where
the slope of the characteristic curve is steepest.

As with other piezoresistive tactile sensors, the choice of pull-up resistance
modulates its response, with a large resistances resulting in a narrow
measuring range and large sensitivity, and vice versa. The ability to adjust
its behavior to capture expected contact intensities for different tasks makes
our design highly versatile.

Similar to other piezoresistive sensors [5, 7, 21], our sensor exhibits
hysteresis, because the EeonTex fabric needs time to relax after contact-based
compression. The larger the hysteresis, the smaller a sensor’s sensibility for
high frequency changes. We measured an average relaxation time of ca.
$0.41\text{\,}\mathrm{s}$ after unloading has ended.

Figure 3 shows measured responses during the loading process for different
pull-up resistors, and the hysteresis during the loading/unloading cycle.
These responses are averaged across ten trials each and highly predictable
with a standard deviation of only $0.29\text{\,}\mathrm{V}$ on average across
contact pressures and $0.67\text{\,}\mathrm{V}$ at maximum. We did not
experience cross-talk at neighboring taxels.

###  IV-B Effect on Compliance

Figure 4: The tactile sensor has only limited effect on finger compliance.
Experimental setup for measuring finger deformations due to horizontal pulling
forces in flexional (top left), extensional (top center), and lateral (top
right) direction with a force-torque sensor which is hooked into the
fingertip. Bottom row: fingertip displacements along horizontal axis for
different pulling forces. Solid lines show mean deformation over 15 trials
(five fingers, three trials per finger), intervals indicate standard
deviations.

A sensor that substantially reduces compliance would impair the grasping and
manipulation capabilities of a soft hand. Therefore, we analyze our sensor’s
effect on compliance when attached to the soft fingers of the RBO Hand 2.

We compare five sensorized to five non-sensorized fingers when applying
pulling forces of $0.5\text{\,}\mathrm{N}$, $1\text{\,}\mathrm{N}$, and
$2\text{\,}\mathrm{N}$ in flexional, extensional, and lateral directions (Fig.
4). To test a finger, it is fixated at its base and a force-torque sensor is
connected to its fingertip via a metal hook. The force-torque sensor is pulled
away from the finger until a specific pulling force is reached while its
vertical position is adjusted to maintain horizontal pulling forces. We
determine the resulting deformations by measuring the distance along the
horizontal axis between the fingertip position in its relaxed and deformed
pose. For each finger, we repeat this process three times for each pulling
force in each direction.

Figure 4 shows the resulting finger deformations for the different forces and
directions. The average distance between fingertip positions is
$11.9\text{\,}\mathrm{m}\mathrm{m}$ for flexional,
$7.7\text{\,}\mathrm{m}\mathrm{m}$ for extensional, and
$8.4\text{\,}\mathrm{m}\mathrm{m}$ for lateral forces. On average, the
sensorized finger deforms about 19.5% less than its non-sensorized
counterpart. The measurable effect on compliance is not surprising since we
added material to the finger. However, this difference does not impair the
hand’s ability to passively adapt its morphology to the shape of grasped
objects, as we show in the following section. The limited effect on
compliance—despite the non-extensibility of the sensor—stems from the fact
that the finger’s palmar side is reinforced by fabric that itself is non-
extensible so that it bends upon inflation. The ability to quantify the loss
of compliance due to sensorization allows making informed adjustments to the
finger’s geometry or to the choice of materials in order to compensate for
this loss, if needed.

##  V In-Hand Object Recognition

We demonstrate our sensor’s ability to provide useful contact information by
enabling a sensorized hand to reliably recognize grasped objects. High
classification rates indicate a sensor’s ability to capture distinctive object
features from tactile information. We compare object recognition performance
based on only tactile sensor information to a baseline case when only air
pressure information from the pneumatic actuators is considered. Furthermore,
we analyze failure cases with respect to physical similarities between objects
and investigate how different levels of spatial resolution contribute to
successful object recognition.

Figure 5: Compliant, sensorized RBO Hand 2 passively adapts its shape to the
grasped object. Left: non-inflated hand. Center: closed hand in absence of an
object after performing the inflation pattern used for all grasps. Right:
hand-object configuration after compliantly grasping a sphere (not part of the
object set).

###  V-A Experimental Setup

We equip the soft pneumatic RBO Hand 2 with tactile sensing by attaching our
sensor to the object-facing (palmar) side of its four silicone-based fingers
(Fig. 5). For comparing our sensor to a baseline case, we also record the
pressure inside each air chamber with Freescale MPX4250 sensors. Each finger
has two of these air chambers: one at its base (proximal) and one at its tip
(distal). The actuators for thumb movement do not provide any data, but are
used for grasping.

We adjust the sensor’s measuring range and sensitivity to match expected
contact intensities. We found that during grasping, contact pressures remain
in more than 97.9% of the cases below 75 kPA. Based on the sensor’s response
behavior (Fig. 3), we choose a pull-up resistance of
$470\text{\,}\mathrm{k}\Omega$ so that these pressures are well within the
sensor’s measuring range.

To evaluate the ability of the sensorized hand to reliably recognize every-day
objects, we chose objects with diverse physical properties that fall into
three groups regarding shape and softness (Fig. 6). Objects within the same
group exhibit highly similar physical properties and are therefore difficult
to distinguish. The three object groups are i) objects that can move inside
deformable packaging, including cat treats, safety glasses, and chopsticks,
ii) rigid plastic bottles, including a whiteboard cleaner, a gum container,
and wet glue, and iii) spherical objects, including a tennis ball, a soft
ball, and a gel ball which all have similar diameters, but different degrees
of softness.

###  V-B Data Acquisition

Figure 6: Object set with nine objects from three object groups, from left to
right: objects in packaging (safety glasses, cat treats, chopsticks), plastic
bottles (whiteboard cleaner, gum container, wet glue), and spherical objects
(soft ball, tennis ball, gel ball).

To record sensory feedback, one of the objects is placed inside the open hand
with its weight resting on the palm which is supported by rigid material,
preventing the object from falling out and its weight from influencing contact
measurements (Fig. 5). Objects are placed in random orientation, resulting in
a variety of contact patterns for the same object. Hand closure brings the
tactile sensors in contact with the object while all of the fingers and also
the thumb are actuated (Fig. 5). The same closing-synergy, based on pre-
defined air masses, is executed for each object. Because of the hand’s
intrinsic compliance, it passively adapts to the shape of the object,
resulting in distinct contact patterns. Once the hand is fully closed, a delay
of three seconds ensures that all sensors have stabilized after inflation-
dependent disturbances. Responses of the tactile sensor with 14 taxels and of
the two air pressure sensors is measured for each of the four fingers,
resulting in a data point of $4\times(14+2)=64$ dimensions. We repeat this
process 30 times for each of the nine objects, resulting in 270 data points in
total.

###  V-C Classification Results

Figure 7: Confusion matrices for object recognition based on tactile (violet)
and air pressure (blue) information: entries indicate the frequency of a label
being assigned to a specific object by the classifier. Results were obtained
by leave-one-out cross-validation. The tactile sensor exhibits superior
performance for all objects except cat treats and chopsticks.

We analyze recorded contacts for enabling the sensorized hand to correctly
classify grasped objects. Specifically, we use machine learning to train a
classifier based on the acquired data set. Training was performed based on
either only tactile information or only air pressure information (56 and 8
dimensions, respectively). To optimize classification accuracy, we compare
different machine learning algorithms and when applicable, perform grid search
for hyperparameter optimization. Tested algorithms include k-nearest neighbors
(KNN), decision trees (DT), and support vector machines (SVM) with linear,
polynomial and radial basis function (RBF) kernels. Classification performance
was evaluated by leave-one-out cross-validation (Table I).

| KNN | DT | SVM (Linear) | SVM (Polynomial) | SVM (RBF)
---|---|---|---|---|---
Tactile | 77.8 | 73.3 | 88.2 | 87.4 | 89.3
Air Pressure | 55.2 | 66.3 | 69.3 | 71.5 | 73.0

TABLE I: Object recognition performance [%] of different classification
algorithms: the SVM with a RBF kernel performs best for both sensor
modalities.

We found that an SVM with an RBF kernel performs best for both sensory
modalities. In the following, we will refer to this classifier when discussing
results (Fig.7). On average, tactile-based classification achieves an accuracy
of 89.3%, which is comparable with performances of other sensorized soft hands
[23, 24]. In comparison, pressure-based classification achieves only 73.0%,
highlighting the superiority of our tactile sensor. However, for cat treats
and chopsticks, air pressure sensing outperforms tactile sensing by 10%, with
83% and 87% compared to 73% and 77%, respectively. Both are objects in
packaging, indicating that classification performance varies across object
groups.

###  V-D Failure Case Analysis

Since classification performance varies across object groups, we now
investigate failure cases in more detail with respect to physical similarities
between objects. For this, we categorize assigned labels into three
categories: i) label is correct, ii) label is incorrect but part of the same
object group as the grasped object, iii) label and object group are incorrect.

The frequencies of these three categories are depicted in Figure 8. On
average, of all incorrectly assigned tactile-based labels, 79.3% are within
the correct object group, compared to only 42.5% of the incorrectly assigned
air pressure-based labels. Thus, not only is tactile sensing superior in
correctly assigning labels, but also in its incorrect labels which are more
often part of the correct object group. In particular, for objects in
packages, for which air pressure-based classification yields slightly better
results, the tactile-based classifier assigns nearly all of the incorrect
labels to the correct object group with 95.3% of the cases, compared to only
36.8% for air pressure sensing.

We argue that the tactile sensor’s performance for objects in packages is
inferior, because these objects have moving parts inside, resulting in a large
variety of tactile responses when grasping the same object multiple times. On
the other side, the coarse sensing of air pressure sensors is less affected by
these movements, as long as object volume and deformation of the compliant
fingers remain similar.

Figure 8: Classification results for tactile (violet) and air pressure (blue)
sensing: assigned labels are either correct (dark color), incorrect but within
the same object group as the grasped object (medium color), or outside this
group (light color). Incorrect tactile-based labels are more frequently within
the correct object group.

###  V-E Spatial Resolution

Figure 9: Tactile sensor with different spatial resolutions, simulated by
averaging over contact information from taxels of the same color, number of
taxels from left to right: one, two, four, eight with vertically combined
fingertip taxels, eight with horizontally combined fingertip taxels, fourteen.

We demonstrate that object classification performance improves with increasing
spatial resolution. For this, we simulate reduced spatial resolutions of
tactile and air pressure sensors by averaging over feedback from multiple
taxels, and from the two air chambers, respectively. The reduced tactile
resolutions result in taxel configurations as depicted in Figure 9, ranging
from a single taxel to fourteen taxels per finger. For the air pressure
sensor, this leads to configurations with a single or with two air chambers.
We trained classifiers based on the different simulated taxel combination,
after optimizing hyperparameters separately for each case.

Object recognition performance for the different spatial resolutions of the
two sensor modalities is shown in Figure 10. Interestingly, the two sensor
modalities achieve comparable performances for the same number of sensing
units. For a single taxel, tactile sensing achieves 55.2% and air pressure
sensing 55.9%. For two taxels, performance rises to 71.5% and 72.9%,
respectively. This indicates the ability of pressure sensing to capture
relevant contact information. However, the number of sensing units of these
sensors is determined by the number of air chambers and their response is
influenced by the inflation state of the soft actuator.

The results indicate that for both sensory modalities, more tactile units
result in better classification performances. This, of course, is not
surprising. However, it demonstrates our sensor’s ability to provide rich
sensory feedback despite its simple design. Furthermore, the results indicate
that also different structure in combining the same number of taxels yields
different results (compare vertically vs. horizontally combined fingertip
taxels in Figure. 10) which is an interesting topic for future research.

Figure 10: Object recognition performance improves for increasing spatial
resolution of each finger’s tactile sensing (violet) and air pressure sensing
(blue). Combining taxel values at the fingertip horizontally (h) yields better
results than vertically (v).

##  VI Conclusion

We proposed a compliant sensor for soft robotics. Its design, based on thin
and flexible materials, minimizes detrimental effects on morphology and
compliance when attached to a soft actuator. The sensor features 14 taxels,
realized by electrodes on a flexible printed circuit board. The sensor can be
manufactured at low-cost in very easy steps within minutes. We characterized
the sensor and explained how choosing an appropriate pull-up resistance leads
to adjustments in sensitivity and measuring range to match task requirements.
The sensor was evaluated on the soft pneumatic RBO Hand 2. In this context, we
showed that the sensor has only little effect on compliance. We also showed
that our sensor provides useful tactile information by demonstrating the
sensorized hand’s capability to identify grasped objects. Finally, we compared
classification performance based on tactile and air pressure information and
showed that tactile sensing recognizes objects more successfully, thanks to
its superior spatial resolution.

## References

  * [1] A. Sieverling, C. Eppner, F. Wolff, and O. Brock, “Interleaving Motion in Contact and in Free Space for Planning Under Uncertainty,” in IEEE/RSJ Int. Conf. on Intelligent Robots and Systems (IROS), pp. 4011–4017, 2017.
  * [2] A. Gupta, C. Eppner, S. Levine, and P. Abbeel, “Learning Dexterous Manipulation for a Soft Robotic Hand from Human Demonstrations,” in IEEE/RSJ Int. Conf. on Intelligent Robots and Systems (IROS), pp. 3786–3793, 2016.
  * [3] E. Páll, A. Sieverling, and O. Brock, “Contingent Contact-Based Motion Planning,” in IEEE/RSJ Int. Conf. on Intelligent Robots and Systems (IROS), pp. 6615–6621, 2018.
  * [4] R. Deimel and O. Brock, “A Novel Type of Compliant and Underactuated Robotic Hand for Dexterous Grasping,” The Int. J. of Robotics Research, vol. 35, no. 1-3, pp. 161–185, 2016.
  * [5] H. Yousef, M. Boukallel, and K. Althoefer, “Tactile sensing for dexterous in-hand manipulation in robotics - A review,” Sensors and Actuators, A: Physical, vol. 167, no. 2, pp. 171–187, 2011.
  * [6] N. Lu and D.-H. Kim, “Flexible and Stretchable Electronics Paving the Way for Soft Robotics,” Soft Robotics, vol. 1, no. 1, pp. 53–62, 2014.
  * [7] Z. Kappassov, J.-A. Corrales, and V. Perdereau, “Tactile sensing in dexterous robot hands — review,” Robotics and Autonomous Systems, vol. 74, pp. 195 – 220, 2015.
  * [8] Y. L. Park, B. R. Chen, and R. J. Wood, “Design and Fabrication of Soft Artificial Skin Using Embedded Microchannels and Liquid Conductors,” IEEE Sensors J., vol. 12, no. 8, pp. 2711–2718, 2012.
  * [9] X. Fang, Z. Liu, Y. Hao, H. Yang, J. Liu, Z. Xie, and L. Wen, “A Soft Actuator with Tunable Mechanical Configurations for Object Grasping Based on Sensory Feedback,” 2019 2nd IEEE Int. Conf. on Soft Robotics (RoboSoft), pp. 25–30, 2019.
  * [10] V. Wall, G. Zöller, and O. Brock, “A Method for Sensorizing Soft Actuators and Its Application to the RBO Hand 2,” in IEEE Int. Conf. on Robotics and Automation (ICRA), pp. 4965–4970, 2017.
  * [11] V. Wall and O. Brock, “Multi-Task Sensorization of Soft Actuators Using Prior Knowledge,” in IEEE Int. Conf. on Robotics and Automation (ICRA), pp. 9416–9421, 2019.
  * [12] G. Zöller, V. Wall, and O. Brock, “Acoustic Sensing for Soft Pneumatic Actuators,” in IEEE/RSJ Int. Conf. on Intelligent Robots and Systems (IROS), pp. 6986–6991, 2018.
  * [13] G. Zöller, V. Wall, and O. Brock, “Active Acoustic Contact Sensing for Soft Pneumatic Actuators,” in IEEE Int. Conf. on Robotics and Automation (ICRA), pp. 7966–7972, 2020.
  * [14] B. Ward-Cherrier, N. Pestell, L. Cramphorn, B. Winstone, M. E. Giannaccini, J. Rossiter, and N. F. Lepora, “The TacTip Family: Soft Optical Tactile Sensors with 3D-Printed Biomimetic Morphologies,” Soft Robotics, vol. 5, no. 2, pp. 216–227, 2018.
  * [15] B. Ward-Cherrier, J. Conradt, M. G. Catalano, M. Bianchi, and N. F. Lepora, “A Miniaturised Neuromorphic Tactile Sensor integrated with an Anthropomorphic Robot Hand,” in 2020 IEEE/RSJ Int. Conf. on Intelligent Robots and Systems (IROS), pp. 9883–9889, 2020.
  * [16] W. Yuan, S. Dong, and E. H. Adelson, “GelSight: High-Resolution Robot Tactile Sensors for Estimating Geometry and Force,” Sensors, vol. 17, no. 12, p. 2762, 2017.
  * [17] M. Lacasse, V. Duchaine, and C. Gosselin, “Characterization of the Electrical Resistance of Carbon-Black-Filled Silicone: Application to a Flexible and Stretchable Robot Skin,” in IEEE Int. Conf. on Robotics and Automation (ICRA), pp. 4842–4848, 2010.
  * [18] S. Stassi, V. Cauda, G. Canavese, and C. F. Pirri, “Flexible Tactile Sensing Based on Piezoresistive Composites: A Review,” Sensors, vol. 14, no. 3, pp. 5296–5332, 2014.
  * [19] J. H. Low, W. W. Lee, P. M. Khin, S. L. Kukreja, H. L. Ren, N. V. Thakor, and C. H. Yeow, “A Compliant Modular Robotic Hand with Fabric Force Sensor for Multiple Versatile Grasping Modes,” in IEEE Int. Conf. on Biomedical Robotics and Biomechatronics (BioRob), pp. 1230–1235, 2016\.
  * [20] J. M. Butt, H. Wang, and R. Pathan, “Design, Fabrication, and Analysis of a Sensorized Soft Robotic Gripper,” 8th Annual IEEE Int. Conf. on Cyber Technology in Automation, Control and Intelligent Systems, CYBER 2018, pp. 169–174, 2019.
  * [21] G. H. Büscher, R. Kõiva, C. Schürmann, R. Haschke, and H. J. Ritter, “Flexible and stretchable fabric-based tactile sensor,” Robotics and Autonomous Systems, vol. 63, no. P3, pp. 244–252, 2015.
  * [22] D. Chen, Y. Cai, and M. C. Huang, “Customizable Pressure Sensor Array: Design and Evaluation,” IEEE Sensors Journal, vol. 18, no. 15, pp. 6337–6344, 2018.
  * [23] C. Jiao, B. Lian, Z. Wang, Y. Song, and T. Sun, “Visual-tactile object recognition of a soft gripper based on faster region-based convolutional neural network and machining learning algorithm,” International Journal of Advanced Robotic Systems, vol. 17, no. 5, p. 1729881420948727, 2020.
  * [24] B. S. Homberg, R. K. Katzschmann, M. R. Dogar, and D. Rus, “Haptic Identification of Objects using a Modular Soft Robotic Gripper,” in IEEE/RSJ Int. Conf. on Intelligent Robots and Systems (IROS), pp. 1698–1705, 2015.


11institutetext: Graduate Institute of Astronomy, National Central University,
Zhongli Dist., Taoyuan City 32001, Taiwan, (R.O.C.)
22institutetext: Department of physics, National Central University, Zhongli
Dist., Taoyuan City 32001, Taiwan, (R.O.C.)
33institutetext: I. E. Tamm Theoretical Physics Division of P. N. Lebedev
Institute of Physics, Leninskii pr. 53, 119991 Moscow, Russia
44institutetext: Moscow Institute of Physics and Technology (State
University), 9, Institutsky lane, Dolgoprudny, 141707, Russia
44email<EMAIL_ADDRESS>44email<EMAIL_ADDRESS>
# Outflows in the presence of cosmic rays and waves with cooling

C. M. Ko 1122 B. Ramzan 11 D. O. Chernyshov 3344

###### Abstract

Context. Plasma outflow from a gravitational potential well with cosmic rays
and self-excited Alfvén waves with cooling and wave damping is studied in the
hydrodynamics regime.

Aims. We study outflows in the presence of cosmic ray and Alfvén waves
including the effect of cooling and wave damping. We seek physically allowable
steady-state subsonic-supersonic transonic solutions.

Methods. We adopted a multi-fluid hydrodynamical model for the cosmic ray
plasma system. Thermal plasma, cosmic rays, and self-excited Alfvén waves are
treated as fluids. Interactions such as cosmic-ray streaming instability,
cooling, and wave damping were fully taken into account. We considered one-
dimensional geometry and explored steady-state solutions. The model is reduced
to a set of ordinary differential equations, which we solved for subsonic-
supersonic transonic solutions with given boundary conditions at the base of
the gravitational potential well.

Results. We find that physically allowable subsonic-supersonic transonic
solutions exist for a wide range of parameters. We studied the three-fluid
system (considering only forward-propagating Alfvén waves) in detail. We
examined the cases with and without cosmic ray diffusion separately.
Comparisons of solutions with and without cooling and with and without wave
damping for the same set of boundary conditions (on density, pressures of
thermal gas, cosmic rays and waves) are presented. We also present the
interesting case of a four-fluid system (both forward- and backward-
propagating Alfvén waves are included), highlighting the intriguing relation
between different components.

###### Key Words.:

Outflows – Hydrodynamics– Cosmic rays – Alfvén waves– Cooling – Wave Damping

##  1 Introduction

Cosmic ray energy density and the energy density of other components of the
interstellar medium (e.g., different phases of gas, magnetic field) are on the
same order of magnitude (e.g., Ferrière, 2001; Cox, 2005). The equipartition
of energy between the different components indicates there are significant
interactions among these components of the interstellar medium. Cosmic rays
can play a dynamical role in the structure and evolution of the interstellar
as well as intergalactic medium. For instances, from the hydrodynamic
perspective, cosmic rays can affect instabilities, such as, Parker
instability, Jeans instability, magnetorotional instability (e.g., Parker,
1969, 1966; Hanasz & Lesch, 2000, 2003; Ryu et al., 2003; Kuwabara et al.,
2004; Kuwabara & Ko, 2006; Ko & Lo, 2009; Lo et al., 2011; Kuwabara & Ko,
2015; Heintz & Zweibel, 2018; Heintz et al., 2020; Kuwabara & Ko., 2020), and
they can modify structures and outflows (e.g., Ko et al., 1991; Yang et al.,
2012; Girichidis et al., 2016; Dorfi & Breitschwerdt, 2012; Recchia et al.,
2016; Ruszkowski et al., 2017; Mao & Ostriker, 2018; Farber et al., 2018;
Holguin et al., 2019; Dorfi et al., 2019; Yu et al., 2020; Recchia, 2020;
Ramzan et al., 2020).

The cosmic-ray driven wind mechanism was first discussed by Ipavich (1975).
Subsequently, Breitschwerdt et al. (1991, 1993), and Zirakashvili et al.
(1996), laid down the proper formulation of cosmic ray propagation with
diffusion and wave damping for cosmic-ray driven winds. Later studies (e.g.,
Everett et al., 2008; Recchia et al., 2016; Wiener et al., 2017; Ruszkowski et
al., 2017; Farber et al., 2018) showed that the continuous gas outflows from
galaxies can affect their evolution directly. For more details, we refer to
the works of Zhang (2018) and Recchia (2020).

Steady-state thermal wind transonic solutions with cooling for a wide range of
mass-loading factors, energy-loading factors, galaxy mass, and galaxy radius
was studied by Bustard et al. (2016). These authors showed that radiative
losses or cooling can enormously affect the wind solutions. The cooling of the
gas is important for finding smooth transition from the subsonic to supersonic
branch (transonic solutions). Bustard et al. (2016) presented the cases of
thermal driven wind with cooling and provided feasible transonic solutions. In
this article, we investigate the effect of cooling on outflows under the
dynamical effect of cosmic rays and waves. We explore different possibilities
of subsonic-supersonic transonic solutions when cosmic ray diffusion and wave
damping mechanism are considered. We note that the terms “subsonic” and
“supersonic” in this article do not refer to the propagation speed of a
disturbance in the medium, but rather the speed that is determined by the
critical point condition of the system. We call this speed the “effective
sound speed.”

The paper is organized as follows. In Section 2, we present the four-fluid
model of the cosmic-ray-plasma system in flux-tube formation. The four fluids
are thermal plasma, cosmic ray, forward- and backward-propagating Alfvén
waves. A simpler system, namely, a system with only a forward-propagating wave
(i.e., a three-fluid system or a one-wave system) is studied in Section 4 and
outflows with and without cosmic ray diffusion are discussed in detail.
Section 5 describes an example of the four-fluid system. In Section 6, we
present our remarks and a summary.

##  2 Hydrodynamical model for a cosmic ray-plasma system

Cosmic rays interact with magnetised plasma via the embedded magnetic field,
in particular, magnetic irregularities or waves. As a result, cosmic rays
advect with and get diffused through the plasma. As cosmic rays stream through
the plasma, they excite hydromagnetic waves via streaming instabilities. When
the energy density of cosmic rays (and waves) are on the same order as that of
the thermal plasma, their feedback onto the thermal plasma cannot be
neglected. Energy exchange is facilitated by various pressure gradients and
wave damping. To study the dynamics of the cosmic ray plasma system, we
adopted the four-fluid hydrodynamical model (or two-wave system) put forward
by (Ko, 1992) (see also Zweibel 2017). In this model, the thermal plasma,
cosmic rays, and two Alfvén waves (one propagating forward along the magnetic
field and one propagating backward) are treated as fluids. In the flux-tube
formulation (i.e., considering the dynamics along the magnetic field only) and
for the steady state, the governing equations are (Ramzan et al., 2020):

$B\Delta=\psi_{B}\,,$ (1) $\rho U\Delta=\psi_{m}\,,$ (2) $\rho
U\frac{dU}{d\xi}=-\frac{d}{d\xi}\left(P_{g}+P_{c}+P^{+}_{w}+P^{-}_{w}\right)+\rho
g_{e}\,,$ (3)
$\frac{1}{\Delta}\frac{dF_{g}\Delta}{d\xi}=U\frac{dP_{g}}{d\xi}+L^{+}_{w}+L^{-}_{w}-\Gamma\,,$
(4)
$\frac{1}{\Delta}\frac{dF_{c}\Delta}{d\xi}=\left[U+(e_{+}-e_{-})V_{A}\right]\frac{dP_{c}}{d\xi}+\frac{P_{c}}{\tau}\,,$
(5)
$\frac{1}{\Delta}\frac{dF^{\pm}_{w}\Delta}{d\xi}=U\frac{dP^{\pm}_{w}}{d\xi}\mp
e_{\pm}V_{A}\frac{dP_{c}}{d\xi}-\frac{P_{c}}{2\tau}-L^{\pm}_{w}\,,$ (6)

where the energy fluxes are

$F_{g}=\left(E_{g}+P_{g}\right)U=\frac{\gamma_{g}P_{g}}{\left(\gamma_{g}-1\right)}\,U\,,$
(7) $\displaystyle F_{c}$
$\displaystyle=\left(E_{c}+P_{c}\right)\left[U+(e_{+}-e_{-})V_{A}\right]-D_{c}$
$\displaystyle=\frac{\gamma_{c}P_{c}}{\left(\gamma_{c}-1\right)}\left[U+(e_{+}-e_{-})V_{A}\right]-\frac{\kappa}{\left(\gamma_{c}-1\right)}\frac{dP_{c}}{d\xi}\,,$
(8) $F^{\pm}_{w}=E^{\pm}_{w}\left(U\pm
V_{A}\right)+P^{\pm}_{w}U=P^{\pm}_{w}\left(3U\pm 2V_{A}\right)\,.$ (9)

Equation (1) and Equation (2) come from the divergent free magnetic field and
mass continuity equation, $\psi_{B}$ and $\psi_{m}$ are called the magnetic
flux and the mass flow rate, respectively; $\Delta(\xi)$ is the cross-
sectional area of the magnetic flux tube and $\xi$ is the coordinate along the
magnetic field. In Equation (3), $\rho$ and $U$ are the density and flow
velocity of the plasma, $P_{g}$, $P_{c}$, and $P^{\pm}_{w}$ are the pressures
of the thermal plasma, cosmic rays, and waves ($\pm$ denote forward- and
backward-propagating Alfvén waves). Also, $g_{e}$ is the external
gravitational field. Equation (4) describes the influence of gas by cooling,
$\Gamma$, and by heating via wave damping, $L^{\pm}_{w}$. The terms $\mp
e_{\pm}V_{A}dP_{c}/d\xi$ in Equations (5) and (6) are wave excitation by
cosmic ray streaming instability. In Equation (2), $\kappa$ is the diffusion
coefficient of cosmic rays and $D_{c}=\kappa dE_{c}/d\xi$ is the cosmic ray
diffusive flux, $P_{c}/\tau$ represents stochastic acceleration. Also,
$V_{A}=B/\sqrt{\mu_{0}\rho}$ is the Alfvén speed. We adopted the following
working model for $e_{\pm}$, $\kappa$ and $1/\tau$ (e.g., Skilling, 1975; Ko,
1992):

$\displaystyle e_{\pm}$ $\displaystyle=\frac{\nu_{\pm}}{(\nu_{+}+\nu_{-})}\,,$
$\displaystyle\kappa$ $\displaystyle=\frac{c^{2}}{3(\nu_{+}+\nu_{-})}\,,$ (10)
$\displaystyle\frac{1}{\tau}$
$\displaystyle=\frac{16\nu_{+}\nu_{-}V_{A}^{2}}{(\nu_{+}+\nu_{-})c^{2}}\,,$

where $\nu_{\pm}$ are the collision frequencies of cosmic rays by forward- and
backward-propagating waves. We take $\nu_{\pm}=c^{2}P^{\pm}_{w}/\eta$ in the
calculations in Sections 3 and 4. The smaller the $\eta$, the stronger the
coupling strength.

##  3 Wind equation in constant flux-tube geometry

In the following sections, we focus on outflows (especially transonic
outflows) under the influence of gas cooling, wave damping, and cosmic ray
diffusion as well. To avoid complications arising from a divergent flux-tube,
we restrict our discussions to constant flux-tube geometry only, that is, we
take $\Delta$ to be a constant. To study outflows, it is useful to express the
momentum equation (Equation (3)) in the form of a wind equation (cf. the
classic stellar wind in Parker (1958)),

$\displaystyle\left(1-M^{-2}_{\rm eff}\right)U\frac{dU}{d\xi}=$
$\displaystyle\quad
g_{e}-\frac{1}{\rho}\frac{dP_{c}}{d\xi}\left[\frac{e_{+}(1+\frac{1}{2}M_{A}^{-1})}{(1+M_{A}^{-1})}+\frac{e_{-}(1-\frac{1}{2}M_{A}^{-1})}{(1-M_{A}^{-1})}\right]$
$\displaystyle\quad+\frac{1}{\rho
U}\left[\frac{P_{c}}{2(1-M_{A}^{-2})\tau}-(\gamma_{g}-1)\left(L^{+}_{w}+L^{-}_{w}-\Gamma\right)\right.$
$\displaystyle\quad\quad\quad\left.+\frac{L^{+}_{w}}{2(1+M_{A}^{-1})}+\frac{L^{-}_{w}}{2(1-M_{A}^{-1})}\right]\,.$
(11)

where
$M_{A}=\sqrt{\mu_{0}}\psi_{m}/(\psi_{B}\sqrt{\rho})=1/({\tilde{\psi}\sqrt{\rho}})$
is the Alfvén Mach number and $M_{\rm eff}=U/a_{\rm eff}$ is the “effective
Mach number.” $dP_{c}/d\xi$ appears explicitly only if cosmic ray diffusion
coefficient is non-zero. The case without cosmic ray diffusion is discussed in
Section 4.1. Here the “effective sound speed” is:

$a^{2}_{\rm
eff}=a^{2}_{g}+a^{+\,2}_{w}\frac{\left(1+\frac{1}{3}M_{A}^{-1}\right)}{\left(1+M_{A}^{-1}\right)}+a^{-\,2}_{w}\frac{\left(1-\frac{1}{3}M_{A}^{-1}\right)}{\left(1-M_{A}^{-1}\right)}\,,$
(12)

where

$a^{2}_{g}=\frac{\gamma_{g}P_{g}}{\rho}\,,\quad
a^{+\,2}_{w}=\frac{3P^{+}_{w}}{2\rho}\,,\quad
a^{-\,2}_{w}=\frac{3P^{-}_{w}}{2\rho}\,.$ (13)

We note that the “effective sound speed” is not necessarily the speed of the
propagation of a disturbance in the medium. It is determined by the “critical
point” of the system, as seen in the term on the left-hand side of Equation
(3); more details are given in Appendix A.

Equations for the gas, cosmic rays, and waves are:

$\frac{dP_{g}}{d\xi}=-\frac{\gamma_{g}P_{g}}{U}\frac{dU}{d\xi}+\frac{(\gamma_{g}-1)}{U}\left(L^{+}_{w}+L^{-}_{w}-\Gamma\right)\,,$
(14)

$\displaystyle\frac{dP^{\pm}_{w}}{d\xi}=-\frac{3P^{\pm}_{w}}{2U}\frac{(1\pm\frac{1}{3}M_{A}^{-1})}{(1\pm
M_{A}^{-1})}\frac{dU}{d\xi}\mp\frac{e_{\pm}M_{A}^{-1}}{2(1\pm
M_{A}^{-1})}\frac{dP_{c}}{d\xi}$ $\displaystyle\quad-\frac{P_{c}}{4U(1\pm
M_{A}^{-1})\tau}-\frac{L^{\pm}_{w}}{2U(1\pm M_{A}^{-1})}\,,$ (15)
$\displaystyle\frac{dP_{c}}{d\xi}=-\frac{\gamma_{c}P_{c}}{U}\frac{[1+\frac{1}{2}(e_{+}-e_{-})M_{A}^{-1}]}{[1+(e_{+}-e_{-})M_{A}^{-1}]}\frac{dU}{d\xi}$
$\displaystyle\quad\quad\quad-\frac{\gamma_{c}P_{c}M_{A}^{-1}}{[1+(e_{+}-e_{-})M_{A}^{-1}]}*\frac{d}{d\xi}(e_{+}-e_{-})$
$\displaystyle\quad\quad\quad+\frac{(\gamma_{c}-1)P_{c}}{U[1+(e_{+}-e_{-})M_{A}^{-1}]\tau}$
$\displaystyle\quad\quad\quad+\frac{1}{U[1+(e_{+}-e_{-})M_{A}^{-1}]}\frac{d}{d\xi}\left(\kappa\frac{dP_{c}}{d\xi}\right)\,.$
(16)

In the following discussion, we take the galactic wind as our fiducial example
of outflows. As a working model, we adopted non-linear Landau damping as the
wave damping mechanism, $L_{w}^{\pm}\propto P_{w}^{\pm 2}\sqrt{P_{g}/\rho
B^{2}}$, and bremsstrahlung radiation for cooling mechanism,
$\Gamma\propto\rho^{2}\sqrt{P_{g}/\rho}$. The corresponding proportionality
constants are denoted by $\delta^{\pm}$ and $\delta_{c}$ in the following
calculations. The gravitational field is $g_{e}=-GM/(\xi+a)^{2}$ (cf. a
Kuzmin-like field near the center of the disk), where $a$ is the
characteristic scale of the field.

##  4 Transonic outflows in a three-fluid system

In this section, we consider a simpler system, namely, a three-fluid or one-
wave system. In this system, only the forward-propagating wave is considered,
namely, $e_{+}=1$, $e_{-}=0$, $1/\tau=0$. We seek physically allowable
solutions from a given set of boundary conditions at the base of the potential
well. Physically allowable solution means the outflow quantities are single-
value with respect to $\xi$, and the velocity and pressures are non-negative
and finite. In particular, we are interested in subsonic-supersonic transonic
solutions. Here, we discuss two different cases: (1) systems without cosmic
ray diffusion and (2) systems with cosmic-ray diffusion. For subsonic or
supersonic solutions, these two cases are not expected be too different
(Ramzan et al., 2020). However, for transonic solutions, they can be very
different because the effective sound speeds for the two cases are very
different (cf. Equation (12) and (18)).

###  4.1 Three-fluid outflow without cosmic ray diffusion

In this subsection, we study three-fluid outflows by neglecting cosmic ray
diffusion (i.e., cosmic ray is strongly coupled to thermal gas). Setting
$e_{+}=1$, $e_{-}=0$, $1/\tau=0$ and $\kappa=0$ in Equation (3), the cosmic-
ray pressure can be expressed analytically as Equation (23) (or see Ramzan et
al. 2020) and the three-fluid wind equation becomes (cf. Equation (3)):

$\displaystyle\left(1-M^{-2}_{\rm
eff}\right)U\frac{dU}{d\xi}=g_{e}+(\gamma_{g}-1)\frac{\Gamma}{\rho U}$ (17)
$\displaystyle+\left[\frac{1}{2(1+M_{A}^{-1})}-(\gamma_{g}-1)\right]\frac{L_{w}^{+}}{\rho
U}\,,$

where the “effective sound speed” is modified to (cf. Equation (12)):

$a^{2}_{\rm
eff}=a^{2}_{g}+a^{+\,2}_{w}\frac{\left(1+\frac{1}{3}M_{A}^{-1}\right)}{\left(1+M_{A}^{-1}\right)}+a^{\,2}_{c}\frac{\left(1+\frac{1}{2}M_{A}^{-1}\right)^{2}}{\left(1+M_{A}^{-1}\right)^{2}}\,,$
(18)

and $a^{2}_{c}=\gamma_{c}P_{c}/\rho$. We explore outflow against gravity,
thus, $g_{e}$ is negative (and approaches zero at large distances away from
the gravitational potential well). For convenience, we call the right hand
side of the wind equation (for this subsection it is Equation (17)) the
driving term. If there is no cooling and no wave damping, then the driving
term is always negative. Thus, in the subsonic regime ($M_{\rm eff}<1$),
$dU/d\xi>0$ and the outflow speed is monotonically increasing with $\xi$. We
may call this accelerating outflow. There are two possibilities: (1) the
accelerating outflow is always subsonic and approaches an asymptotic value at
large distances or (2) the flow speed of the solution approaches the effective
sound speed at finite distance, then $dU/d\xi\rightarrow\infty$ and the
solution becomes unphysical. By the same token, in the supersonic regime
($M_{\rm eff}>1$), $dU/d\xi<0$, and the outflow is a decelerating outflow.
Here, there are also two possibilities: (1) a decelerating supersonic outflow
that approaches a finite velocity at large distances or (2) the solution
approaches the effective sound speed and becomes unphysical.

If cooling is included (but still ignoring wave damping), then the driving
term can be positive or negative. Thus, a new possibility arises. It is then
possible to have a solution that as its flow speed approaches the effective
sound speed and driving term also approaches zero simultaneously. $dU/d\xi$
can be finite at this point (the sonic point or critical point), and the
solution can pass through the critical point. A transonic solution is born. A
more detailed analysis of this possibility can be found in Appendix A.1.

When wave damping is included, its contribution to the driving term can be
positive or negative. For sub-Alfvénic flow the contribution is always
negative if $\gamma_{g}>5/4$, while for super-Alfvénic flow the contribution
is always negative if $\gamma_{g}>3/2$ and always positive if
$\gamma_{g}<5/4$. In this work we only consider cases with $\gamma_{g}=5/3$.
We seek physically allowable solutions to the wind equation and pressure
equations (Equations (14)-(17)) subjected to a set of boundary conditions of
$(U,P_{g},P_{c},P_{w}^{+},\rho,B)$ at the base of the gravitational potential
well. However, not every set of boundary conditions gives us a physically
allowable solution. Physically allowable solutions are divided into subsonic
solutions, supersonic solutions, and transonic solutions. For the purposes of
this study, we are interested in subsonic-supersonic transonic solutions only.
At the outset we chose a set $(P_{gb},P_{cb},P_{wb}^{+},\rho_{b},B_{b})$ at
the base of the potential well and then adjusted $U_{b}$ to obtain the
transonic solution.

In Figure 1, we show a typical result of three-fluid outflows with cooling. We
present different cases side by side. The middle column is the case without
wave damping and the right column is the case with wave damping. The boundary
conditions at the base $(P_{gb},P_{cb},P_{wb}^{+},\rho_{b},B_{b})$ for both
cases are the same, but $U_{b}$ are different in order to obtain a transonic
solution for the two cases. As shown in the figure, the velocity and pressure
profiles of the two cases are quite similar. We also included a pure thermal
outflow for comparison (the left column). The boundary condition for the pure
thermal outflow has the same $(P_{gb},\rho_{b},B_{b})$ as the three-fluid
cases, while $U_{b}$ is adjusted to get the transonic solution. The velocity
of the pure thermal outflow is significantly smaller than the three-fluid
cases (black solid line in the upper row of Figure 1 is flow speed). This is
expected as the corresponding effective sound speeds are different (thin brown
line in the upper row of Figure 1 is effective sound speed). The effective
sound speed of the three-fluid cases has extra contributions from waves and
cosmic rays (see Equation (18)). Nevertheless, the thermal gas pressure
profiles of the three cases are shown to be rather similar (see the lower row
of Figure 1).

Figure 1:  Transonic solutions of a three-fluid system without cosmic ray
diffusion. The transonic solution is subsonic-supersonic transition. In this
figure, $\gamma_{g}=5/3$, $\gamma_{c}=4/3$, $GM=1$, $a=1$. Left column: Pure
thermal outflow with cooling ($\delta_{c}=0.1$). Middle column: Three-fluid
outflow with cooling ($\delta_{c}=0.1$) but without wave damping
($\delta^{+}=0$). Right column: Three-fluid outflow with cooling
($\delta_{c}=0.1$) and wave damping ($\delta^{+}=0.1$). Upper row: Black line
and thin brown lines show the flow speed and the effective sound speed of the
system. For the three-fluid cases, Alfvén speed is also shown (thin gray
lines). Lower row: Blue, orange, and green lines show the pressures of the
thermal gas, cosmic ray, and forward-propagating wave, respectively. The
boundary conditions $(P_{gb},\rho_{b},B_{b})$ are the same for all three
cases, and $(P_{cb},P_{wb}^{+})$ are the same for the two three-fluid cases.
In this example, $P_{gb}=1$, $P_{cb}=0.5$, $P_{wb}^{+}=0.25$, $\rho_{b}=1$,
$B_{b}=\psi_{B}^{\prime}=0.25$. The boundary condition $U_{b}$ is adjusted to
obtain the transonic solutions.

###  4.2 Three-fluid outflow with cosmic ray diffusion

An important feature of cosmic ray transport is the diffusion of cosmic rays
through the thermal plasma. The diffusion coefficient indicates the coupling
of cosmic rays with the plasma; the larger the diffusion coefficient, the
weaker the coupling. We modeled the diffusion coefficient as inversely
proportional to the wave pressure (thus, as the wave pressure tends to zero,
cosmic rays tend to be decoupled from the plasma). The wind equation is
Equation (3) with only the forward-propagating wave, and $e_{+}=1$, $e_{-}=0$,
$1/\tau=0$. We note that cosmic ray pressure gradient appears explicitly in
the driving term of the wind equation (the right hand side of Equation (3)).
In comparison with the cases without cosmic ray diffusion (Section 4.1), the
presence of cosmic ray diffusion shifts the narrative of the existence of
transonic solution in two ways: the cosmic ray contribution to the effective
sound speed disappears (see Equation (12)), and the role of cosmic rays in the
driving term. In fact, transonic solution can exist in three-fluid system with
diffusion even without cooling and wave damping (see Appendix A.2). However,
we have opted to focus on the role of cooling on transonic solution in the
present study.

We explore solutions for the wind equation (Equation (3)), together with the
corresponding pressure equations, namely, Equations (14) to (16) that are
subject to boundary conditions at the base of the gravitational potential
well. In addition to $(U_{b},P_{gb},P_{cb},P_{wb}^{+},\rho_{b},B_{b})$, we
need one more boundary condition for the gradient of cosmic-ray pressure (or
cosmic-ray diffusive flux $D_{cb}$) at the base. In Figure 2, we compare the
case with both cooling and wave damping are “OFF” to cases where wave damping
is “ON” or cooling is “ON,” or both are “ON.” In this example, wave damping is
“OFF” and “ON,” referring to $\delta^{+}=0$ and 0.1 and cooling is “OFF” and
“ON,” referring to $\delta_{c}=0$ and 0.1. All cases have the same boundary
conditions $(U_{b},P_{gb},P_{cb},P_{wb}^{+},\rho_{b},B_{b})$ except $D_{cb}$.
For those for which cooling is “ON,” we adjust $D_{cb}$ to obtain the
transonic solution (subsonic-supersonic transition), while for those for which
cooling is “OFF,” we again adjust $D_{cb}$ to obtain the physically allowable
solution (i.e., subsonic solution that extends to large distances).

In Figure 2, the upper row shows the outflow speed (black), effective sound
speed (thin brown) and Alfvén speed (thin gray). The lower row shows the
thermal pressure (blue), cosmic ray pressure (orange) and wave pressure
(green). The case where both wave damping and cooling are “OFF” is compared to
the case that only wave damping is “ON” in the left column, as well as to the
case that only cooling is “ON” in the middle column, and to the case that both
wave damping and cooling are “ON” in the right column.

In the left column, all the cases are without cooling and they all belong to
the subsonic outflow. The dotted line is the case without wave damping and
cooling. The dashed, dot-dashed, and solid lines belong to wave damping case
and all with the same $\delta^{+}=0.1$, but with slightly different $D_{cb}$.
Here, $D_{cb}$ decreases from dashed to dot-dashed to solid and the solid line
has the minimum $D_{cb}$ beyond which there is no physically allowable
solution in this example. The solid line is barely different from the dotted
line. The dashed and dot-dashed lines have a rapid increase in $U$ and $P_{c}$
(and a rapid decrease in $P_{g}$) somewhere along the flow where $P_{w}^{+}$
becomes diminish and cosmic rays decouple from the thermal gas. These are
typical quasi-thermal outflows discussed in Ramzan et al. (2020). Appendix B
gives a summary of different flow profiles in three-fluid system (with
forward-propagating wave).

In the middle and the right columns, cooling is included and the solutions are
subsonic-supersonic transonic flows (solid lines). The profiles carry the
flavor of quasi-thermal outflows. It appears that the transition (sonic point)
takes place at a further distance for the case without wave damping (the
middle column). However, the profiles at large distances for both cooling
cases are quite similar in this example.

We are also interested in the effect of the cosmic-ray coupling strength on
the results with cooling (in which the physically allowable solutions are
transonic solutions). The cosmic ray diffusion coefficient in the three-fluid
case is $\kappa=\eta/3/P_{w}^{\pm}$. The coupling is stronger for smaller
$\eta$. In Figure 2, $\eta=1$. For the cases with cooling in Figure 2, we work
out the corresponding solutions for $\eta=5$, 2, 1, and 0.5 (by adjusting
$D_{cb}$), and the result are shown in Figure 3. In the figure, the left
column is the case with cooling “ON” and wave damping “OFF” ($\delta_{c}=0.1$,
$\delta^{+}=0$) and the right column with both cooling and wave damping “ON”
($\delta_{c}=0.1$, $\delta^{+}=0.1$). The rapid increase in the flow speed and
cosmic ray pressure (and decrease in the thermal pressure) become steeper when
$\eta$ is smaller. The solutions will be unphysical for even smaller $\eta$.
Moreover, we note that for $\eta=0$ (strong coupling limit), the system
belongs to Section 4.1. The corresponding solution with the same boundary
conditions $(U_{b},P_{gb},P_{cb},P_{wb}^{+},\rho_{b},B_{b})$ as those in
Figure 2 is unphysical.

Figure 2:  Comparison of transonic and subsonic solutions of a three-fluid
system with cosmic ray diffusion. In this figure, $\gamma_{g}=5/3$,
$\gamma_{c}=4/3$, $GM=1$, $a=1$, $\eta=1$. Cases with wave damping and/or
cooling are compared to the cases without both wave damping and cooling. All
cases have the same boundary conditions
$(U_{b},P_{gb},P_{cb},P_{wb}^{+},\rho_{b},B_{b})$ at the base except $D_{cb}$.
The diffusive flux (or cosmic ray pressure gradient) at the base is adjusted
to get physically allowable solutions (subsonic solutions for the cases
without cooling and transonic solution for the cases with cooling). In this
example, $U_{b}=1$, $P_{gb}=2.5$, $P_{cb}=1$, $P_{wb}^{+}=0.15$,
$\rho_{b}=2/3$, $B_{b}=\psi_{B}^{\prime}=1/3$. Dotted line displays the case
in which both wave damping and cooling are “OFF” ($\delta_{c}=0$,
$\delta^{+}=0$). Left column: Dashed, dot-dashed, and solid lines are the
cases in which only wave damping is “ON” ($\delta_{c}=0$, $\delta^{+}=0.1$).
The three solutions correspond to slightly different $D_{cb}$ (see text).
Middle column: Solid line is the case with only cooling is “ON”
($\delta_{c}=0.1$, $\delta^{+}=0$). Right column: Solid line is the case in
which both cooling and wave damping are “ON” ($\delta_{c}=0.1$,
$\delta^{+}=0.1$). Upper row: Black, thin brown, and grey lines show the flow
speed, the effective sound speed and the Alfvén speed, respectively. Lower
row: Blue, orange, and green lines show the pressures of the thermal gas,
cosmic ray, and forward-propagating wave, respectively.

Figure 3:  Comparison of results with different cosmic ray coupling strength.
The cosmic ray diffusion coefficient is defined as
$\kappa=\eta/3/P_{w}^{\pm}$. Here, we compare solutions for different $\eta$
of the cases with cooling in Figure 2. Left column: Cooling is “ON” and wave
damping is “OFF” ($\delta_{c}=0.1$, $\delta^{+}=0$). Right column: Both
cooling and wave damping are “ON” ($\delta_{c}=0.1$, $\delta^{+}=0.1$). Four
different cases are shown, (1) $\eta=5$ (dotted line), (2) $\eta=2$ (dashed
line), (3) $\eta=1$ (dot-dashed line), and (4) $\eta=0.5$ (solid line). All
other parameters are the same as Figure 2, except that $D_{cb}$ is adjusted to
produce transonic solutions.

##  5 Example of a four-fluid outflow

The four-fluid model is the most comprehensive hydrodynamic model for cosmic
ray-plasma system. It involves both forward- and backward-propagating waves.
It is far more complicated than the three-fluid model discussed in Section 4.
When both waves exist, stochastic acceleration cannot be neglected
($1/\tau>0$). In addition, it is possible for the cosmic ray streaming
instability (represented by the cosmic ray pressure gradient term in Equation
(3)) to have an opposite effect on the two waves. Here, we present an
interesting example to illustrate the intriguing relation between different
factors. In Figure 4, we compare the case with cooling to the case without
cooling. Both cases have the same wave damping and the same boundary
conditions $(U_{b},P_{gb},P_{cb},P_{wb}^{+},P_{wb}^{-},\rho_{b},B_{b})$ at the
base of the gravitational potential well. By adjusting the cosmic ray
diffusive flux $D_{cb}$ at the base, we solve the wind equation (Equation (3))
together with the corresponding pressure equations (Equations (14)-(16)) for
transonic solution for the case with cooling and for subsonic solution that
can extend to large distances (i.e., physically allowable solution) for the
case without cooling.

The case with cooling has a number of distinct features, as shown in Figure 4.
The velocity profile (black solid line in the left panel of Figure 4) has an
obvious hump. Non-monotonicity is not observed in cases without cooling. The
rapid decrease of velocity is related to the (second) increase of cosmic ray
pressure (orange solid line in the right panel) and the rise in the pressure
of the backward-propagating wave (purple solid line in the right panel).
Moreover, the backward-propagating wave survives at large distances because
wave damping approaches to zero. Cooling causes the temperature ( $T\propto
P_{g}/\rho$) to decrease. At this point, we recall the working model for
cooling ($\Gamma\propto\rho^{2}\sqrt{P_{g}/\rho}$) and wave damping
($L_{w}^{\pm}\propto P_{w}^{\pm\,2}\sqrt{P_{g}/\rho B^{2}}$). For details, see
the end of Section 2. When the temperature becomes very low, both cooling and
wave damping becomes ineffective.

It is worth pointing out that in this particular example, the wave pressures
of the case with cooling are practically zero between $\xi=10$ and $30$. Thus,
cosmic ray diffusion coefficient becomes very large and cosmic ray is
decoupled from the plasma. However, when the backward-propagating wave rises
again (around $\xi=40$), cosmic rays are recoupled to the plasma again.

Another interesting result of this particular example is the flow speed at
large distances of the subsonic solution is actually larger than the transonic
solution (black dashed and solid lines in the left panel of Figure 4).

Figure 4:  Example of transonic and subsonic solutions of a four-fluid system
with both forward- and backward-propagating waves. In this figure,
$\gamma_{g}=5/3$, $\gamma_{c}=4/3$, $GM=1$, $a=1$, $\eta=1$. Two cases are
compared, (1) dashed line: Cooling is “OFF” ($\delta_{c}=0$) and (2) solid
line: Cooling is “ON” ($\delta_{c}=0.1$). Wave damping is “ON” for both cases
($\delta^{\pm}=1$). Left panel: Black, thin brown, and grey lines show the
flow speed, the effective sound speed, and the Alfvén speed, respectively.
Right panel: Blue, orange, green, and purple lines show the pressures of
thermal gas, cosmic ray, forward-propagating wave, and backward-propagating
wave, respectively. The boundary conditions
$(U_{b},P_{gb},P_{cb},P_{wb}^{+},P_{wb}^{-},\rho_{b},B_{b})$ are the same for
the two cases. In this example, $U_{b}=2/3$, $P_{gb}=1$, $P_{cb}=1$,
$P_{wb}^{+}=0.5$, $P_{wb}^{-}=0.4$, $\rho_{b}=1$,
$B_{b}=\psi_{B}^{\prime}=1/3$. The cosmic ray pressure gradient or diffusive
flux at the base is adjusted to get the transonic solution if cooling is “ON”
and get a physically allowable subsonic solution if the cooling is “OFF.”

Figure 5:  Comparison of solutions of a three-fluid system with cooling but
without cosmic ray diffusion and wave damping. In this figure,
$\gamma_{g}=5/3$, $\gamma_{c}=4/3$, $GM=1$, $a=1$, $\delta_{c}=0.1$,
$\delta^{+}=0$. Four solutions are shown: (1) Solid line: Subsonic-supersonic
transonic solution; (2) Long-dashed line: Supersonic-subsonic transonic
solution; (3) Dot-dashed line: Subsonic solution; and (4) Dashed line:
Supersonic solution. The three-fluid case is shown in the middle and right
columns, with boundary conditions
$(P_{gb},P_{cb},P_{wb}^{+},\rho_{b},B_{b})=(1,0.5,0.25,1,0.25)$. For
comparison, a corresponding pure thermal case is shown in the left column with
boundary conditions $(P_{gb},\rho_{b},B_{b})$=$(1,1,0.25)$. In the top-left
and top-middle panels, the orange and brown dots are the sonic points of the
subsonic-supersonic and supersonic-subsonic transonic solutions, respectively.

##  6 Summary and discussion

In this article, we study the role of cosmic rays and waves on outflows with
cooling. Our fiducial outflow example is that of galactic winds. For
simplicity, we only considered steady-state outflows in a constant flux tube
(i.e., one-dimensional flow). For a pure thermal outflow in subsonic regime,
cooling causes the outflow speed to decrease towards small values and zero at
finite distances. The density becomes very large at finite distances and it
would be difficult to match the density, for instance, in the intergalactic
medium. However, we find that if cosmic rays and waves are present and the
cosmic ray diffusion can be neglected, then still it is possible to have an
outflow to large distances even in the subsonic regime. Let us examine the
case of three-fluid system with no cosmic ray diffusion and no wave damping in
further detail. In this case, $P_{c}$ and $P_{w}^{+}$ can be expressed in
terms of $U$ explicitly (see Equations (23) and (24)). To facilitate
discussion, let us write down the equations governing the flow speed, thermal
pressure, and temperature as:

$\displaystyle\left(1-M^{-2}_{\rm
eff}\right)U\frac{dU}{d\xi}=-\left(g_{e}^{\dagger}-\Gamma^{\dagger}\right)\,,$
(19) $\displaystyle\left(1-M^{-2}_{\rm
eff}\right)\frac{1}{\rho}\frac{dP_{g}}{d\xi}=M^{-2}_{\rm
g}\left(g_{e}^{\dagger}-\beta_{P}\Gamma^{\dagger}\right)\,,$ (20)
$\displaystyle\left(1-M^{-2}_{\rm
eff}\right)\frac{k_{B}}{m}\frac{dT}{d\xi}=\frac{(\gamma_{g}-1)M^{-2}_{\rm
g}}{\gamma_{g}}\left(g_{e}^{\dagger}-\beta_{T}\Gamma^{\dagger}\right)\,,$ (21)

where $M_{\rm eff}=U/a_{\rm eff}$, $a_{\rm eff}$ is given by Equation (18),
$M_{g}^{-2}=\gamma_{g}P_{g}/\rho U^{2}$, $g_{e}^{\dagger}=-g_{e}>0$,
$\Gamma^{\dagger}=(\gamma_{g}-1)\Gamma/\rho U>0$, and

$\displaystyle\beta_{P}=\frac{\left(1-M_{\rm
eff}^{-2}+M^{-2}_{g}\right)}{M^{-2}_{g}},\quad\beta_{T}=\frac{\gamma_{g}\beta_{P}-1}{(\gamma_{g}-1)}\,.$
(22)

We note that $M_{\rm eff}=M_{g}$ for pure thermal outflow. In the supersonic
regime, $\beta_{T}>\beta_{P}>1$, while in the subsonic regime
$\beta_{T}<\beta_{P}<1$ ($\beta_{P}$ can be negative, and $\beta_{T}<0$ if
$\beta_{P}<1/\gamma_{g}$).

In the subsonic regime ($1-M_{\rm eff}^{-2}<0$):

  * •

(a1): If $\beta_{T}\Gamma^{\dagger}$ ¡ $\beta_{P}\Gamma^{\dagger}$ ¡
$\Gamma^{\dagger}$ ¡ $g_{e}^{\dagger}$, then $dU/d\xi$ ¿ 0, $dP_{g}/d\xi$ ¡ 0,
$dT/d\xi$ ¡ 0

  * •

(a2): If $\beta_{T}\Gamma^{\dagger}$ ¡ $\beta_{P}\Gamma^{\dagger}$ ¡
$g_{e}^{\dagger}$ ¡ $\Gamma^{\dagger}$, then $dU/d\xi$ ¡ 0, $dP_{g}/d\xi$ ¡ 0,
$dT/d\xi$ ¡ 0

  * •

(a3): If $\beta_{T}\Gamma^{\dagger}$ ¡ $g_{e}^{\dagger}$ ¡
$\beta_{P}\Gamma^{\dagger}$ ¡ $\Gamma^{\dagger}$, then $dU/d\xi$ ¡ 0,
$dP_{g}/d\xi$ ¿ 0, $dT/d\xi$ ¡ 0

  * •

(a4): If $g_{e}^{\dagger}$ ¡ $\beta_{T}\Gamma^{\dagger}$ ¡
$\beta_{P}\Gamma^{\dagger}$ ¡ $\Gamma^{\dagger}$, then $dU/d\xi$ ¡ 0,
$dP_{g}/d\xi$ ¿ 0, $dT/d\xi$ ¿ 0

In the supersonic regime ($1-M_{\rm eff}^{-2}>0$):

  * •

(b1): If $\Gamma^{\dagger}$ ¡ $\beta_{P}\Gamma^{\dagger}$ ¡
$\beta_{T}\Gamma^{\dagger}$ ¡ $g_{e}^{\dagger}$, then $dU/d\xi$ ¡ 0,
$dP_{g}/d\xi$ ¿ 0, $dT/d\xi$ ¿ 0

  * •

(b2): If $\Gamma^{\dagger}$ ¡ $\beta_{P}\Gamma^{\dagger}$ ¡ $g_{e}^{\dagger}$
¡ $\beta_{T}\Gamma^{\dagger}$, then $dU/d\xi$ ¡ 0, $dP_{g}/d\xi$ ¿ 0,
$dT/d\xi$ ¡ 0

  * •

(b3): If $\Gamma^{\dagger}$ ¡ $g_{e}^{\dagger}$ ¡ $\beta_{P}\Gamma^{\dagger}$
¡ $\beta_{T}\Gamma^{\dagger}$, then $dU/d\xi$ ¡ 0, $dP_{g}/d\xi$ ¡ 0,
$dT/d\xi$ ¡ 0

  * •

(b4): If $g_{e}^{\dagger}$ ¡ $\Gamma^{\dagger}$ ¡ $\beta_{P}\Gamma^{\dagger}$
¡ $\beta_{T}\Gamma^{\dagger}$, then $dU/d\xi$ ¿ 0, $dP_{g}/d\xi$ ¡ 0,
$dT/d\xi$ ¡ 0

We note that a system can exist in a range (a4) of the subsonic regime only
for a small extent of $\xi$ (if at all). In range (a4), as $\xi$ increases,
$U$ decreases, and $P_{g}$ increases, which implies that $\beta_{P}$
decreases. When $\beta_{P}$ is less than $1/\gamma_{g}$, the value of
$\beta_{T}$ is negative and the system will leave the current range (a4) and
enter range (a3).

In the subsonic regime when cooling is present, the flow speed decreases as
$\xi$ increases if the system is in range (a2), (a3) or (a4). To prevent the
flow from stalling at finite distance, the cooling has to be turned off before
this happens. The system will then stay in range (a1), and the flow speed will
gradually increase towards a finite value as $g_{e}\approx 0$ at large
distances. Suppose the cooling function vanishes as temperature goes to a
small value (e.g., the cooling function $\Gamma\propto\rho^{2}\sqrt{T}$
approaches zero as $T$ approaches zero). When the flow is in the range of (a2)
or (a3) of the subsonic regime, it has a tendency to become stalled. In these
ranges, the temperature also decreases. When the temperature decreases rapidly
enough, the cooling can be turned off before the speed of the flow reaches
zero. It is then possible to have a physically allowable solution that extends
to large distances. We then consider what the role of cosmic rays and waves is
on the flow in the subsonic regime. We thus denote $\beta^{\rm therm}_{p}$ and
$\beta^{\rm therm}_{T}$ for the pure thermal case and $\beta^{\rm 3f}_{p}$ and
$\beta^{\rm 3f}_{T}$ for the three-fluid case (see Equation (22) for the
definition of $\beta_{P}$ and $\beta_{T}$). If $U$ and $P_{g}$ are the same in
the two cases, then $\beta_{P}^{\rm therm}>\beta_{P}^{\rm 3f}$ and
$\beta_{T}^{\rm therm}>\beta_{T}^{\rm 3f}$. When $\beta_{T}$ is smaller, the
temperature decreases more rapidly. Hence, the presence of cosmic rays and
waves helps the temperature to drop faster and turns off the cooling. As a
result, it is easier to have a subsonic flow when cosmic rays and waves are
included.

Figure 5 illustrates the effect of cosmic rays and waves on outflows, in
particular, subsonic flows. Radiative cooling by bremsstrahlung
($\Gamma\propto\rho^{2}\sqrt{T}$) is adopted here. The left column shows the
pure thermal case. The middle and right columns show the three-fluid case
(without cosmic ray diffusion and wave damping). The figure shows four types
of solutions, subsonic, subsonic-supersonic, supersonic-subsonic, and
supersonic. All the given solutions have the same boundary conditions
$(P_{gb},\rho_{b},B_{b})$ for both cases and, in addition, the same
$(P_{cb},P_{wb}^{+})$ in the case of three-fluid systems. Also, $U_{b}$ is
adjusted to get the two transonic solutions. $U_{b}$ for the subsonic
(supersonic) solution is slightly smaller (larger) than the one for the
subsonic-supersonic (supersonic-subsonic) transonic solution. As shown in the
top-left panel of the figure for the pure thermal outflow, the subsonic (dot-
dashed line) and the supersonic-subsonic (long dashed line) solutions are not
capable of being extended to large distances and the flow velocity goes to
zero at finite distance. The top-middle panel shows that with the help of
cosmic rays and waves, the temperature drops to small value and the cooling is
off before the flow velocity goes to zero. The flow is then able to enter
range (a1) and approaches a finite flow speed at large distances (see middle
column of Figure 5. The sudden change in $U$ (and $P_{c}$, $P_{w}^{+}$) in the
subsonic region occurs when the temperature drops to zero and the cooling
function also becomes zero. In this particular case, the flow goes from range
(a2) to range (a1).

The orange dot in the top-left and top-middle panels of Figure 4 is the sonic
point of the subsonic-supersonic transonic solution, i.e., the location where
the flow changes from subsonic to supersonic. At this point
$\beta_{T}=\beta_{P}=1$, and the flow moves from range (a1) in the subsonic
regime to range (b4) in the supersonic regime. The brown dot is the sonic
point of the supersonic-subsonic transonic solution, and the flow moves from
range (b1) in the supersonic regime to range (a4) in the subsonic regime (and
leaves range (a4) soon). The apparent crossing of the two transonic solutions
in the figure is only a projection effect.

We conclude that in the case of three-fluid system without cosmic ray
diffusion and wave damping, a subsonic outflow is possible provided that the
amount of cosmic ray and wave is large enough. We caution that if cosmic ray
diffusion is allowed, the line of reasoning above does not apply and the
existence of subsonic flow will be in question (we have tried a range of
parameters for physically allowable subsonic solutions but without success).

As discussed above and shown in Figure 5, a three-fluid subsonic outflow is
possible. However, the flow speed at large distances is small, that is, the
density is large (e.g., the example in Figure 5 has a density at large
distances larger than at the base of the potential well). We are thus more
interested in subsonic-supersonic transonic solutions.

In summary, in this article, we present a study of the steady-state subsonic-
supersonic transonic outflow of cosmic-ray plasma systems with cooling. We
only considered one-dimensional flow. Our fiducial model is based on the
galactic wind. We choose bremsstrahlung radiation as working model for cooling
and nonlinear Landau damping for wave damping.

We examined the three-fluid model or one-wave model (which comprises thermal
plasma, cosmic ray, and self-excited forward-propagating Alfvén wave) without
cosmic ray diffusion in detail (Section 4.1). Here
$(P_{g},P_{c},P_{w}^{+},\rho,B)$ are set at the base of the gravitational
potential well. Here, the value of $U$ at the base is adjusted to obtain a
transonic solution. With the same pressures, density, and magnetic field at
the base, the profiles of the flow speed and pressures of the cases with and
without wave damping are rather similar. However, even with the same thermal
pressure and density at the base, the flow speed of the pure thermal transonic
outflow is significantly lower than the three-fluid flow, which is expected,
as the effective sound speeds are different for both cases (see Figure 1).

When cosmic ray diffusion is taken into consideration in the three-fluid
system, the behavior exhibited by the solutions can be quite different,
particularly when the cosmic-ray diffusion coefficient depends on waves
(Section 4.2). In this case, $(U,P_{g},P_{c},P_{w}^{+},\rho,B)$ are set at the
base and the cosmic ray pressure gradient or diffusive flux ($D_{cb}$) is
adjusted to get transonic solutions. We compare the case of without both
cooling and wave damping to cases of either cooling or wave damping, or both.
In this example, cases with either cooling or wave damping show “rapid change”
in velocity and pressure profiles. This is a characteristic feature of quasi-
thermal outflows (Ramzan et al. 2020; see Appendix B). When cooling is
considered, subsonic-supersonic transition is the only physical solution. In
the presence of wave damping, the “sonic point” seems to be pulled closer to
the base (compare the middle and the right columns of Figure 2). When the
coupling between cosmic rays and the thermal gas increases ($\eta$ decreases),
the “rapid change” in the profiles becomes “steeper” (see Figure 3). In this
example, there will be no physical solution for a small-enough value for
$\eta$.

In Section 5, we gave an example on the most comprehensive cosmic ray-MHD
system, the four-fluid model or two-wave model (where both forward- and
backward-propagating waves are included). Due to complicated interactions
between different components, the profiles of the flow and pressures may show
some interesting features. For instance, in this particular example, when
cooling is “ON” (wave damping is also “ON”) the flow speed has a hump that is
not observed in cases without cooling, Moreover, along the flow cosmic ray is
first coupled to the plasma, subsequently decoupled and then recoupled again
at the further down the flow (indicated by the amount of waves along the flow,
shown as green and purple solid lines in the right panel of Figure 4). In
addition, in this particular example, the flow speed at large distances of the
subsonic solution can be larger than that of the subsonic-supersonic transonic
solution (left panel of Figure 4).

###### Acknowledgements.

CMK and BR are supported in part by the Taiwan Ministry of Science and
Technology grant MOST 109-2112-M-008-005 and MOST 110-2112-M-008-005.

## References

  * Breitschwerdt et al. (1991) Breitschwerdt, D., McKenzie, J. F., & Völk, H. J. 1991, A&A, 245, 79
  * Breitschwerdt et al. (1993) Breitschwerdt, D., McKenzie, J. F., & Völk, H. J. 1993, A&A, 269, 54
  * Bustard et al. (2016) Bustard, C., Zweibel, E. G., & D’Onghia, E. 2016, ApJ, 29, 819
  * Chevalier (1985) Chevalier, R. A.,& Clegg, A. W. 1985, Natur, 317, 44
  * Cox (2005) Cox, D. P. 2005, ARA&A, 43, 1
  * Dorfi & Breitschwerdt (2012) Dorfi, E. A., & Breitschwerdt, D. 2012, A&A, 540, A77
  * Dorfi et al. (2019) Dorfi, E. A., Steiner, D., Ragossnig, F., & Breitschwerdt, D. 2019, A&A, 630, A107
  * Everett et al. (2008) Everett, J. E., Zweibel, E. G., Benjamin, R. A., et al. 2008, ApJ, 674, 258
  * Farber et al. (2018) Farber, R., Ruszkowski, M., Yang, H.-Y. K., & Zweibel, E. G. 2018, ApJ, 856, 112
  * Ferrière (2001) Ferrière, K. M. 2001, Rev. Mod. Phys. 73, 1031
  * Girichidis et al. (2016) Girichidis, P., Naab, T., Walch, S., et al. 2016, ApJ, 816, L19
  * Hanasz & Lesch (2000) Hanasz, M., & Lesch, H. 2000, ApJ, 543, 235
  * Hanasz & Lesch (2003) Hanasz, M., & Lesch, H. 2003, A&A, 412, 331
  * Heintz & Zweibel (2018) Heintz, E., & Zweibel, E. G. 2018, ApJ, 860, 97
  * Heintz et al. (2020) Heintz, E., Bustard, C., & Zweibel, E. G. 2020, ApJ, 891, 157
  * Holguin et al. (2019) Holguin, F., Ruszkowski, M., Lazarian, A., Farber, R., & Yang, H.-Y. K. 2019, MNRAS, 490, 1271
  * Ipavich (1975) Ipavich, F. M. 1975, ApJ, 196, 107
  * Ko (1992) Ko, C. M. 1992, A&A, 259, 377
  * Ko et al. (1991) Ko, C. M., Dougherty, M. K., & McKenzie, J. F. 1991, A&A, 241, 62
  * Ko & Lo (2009) Ko, C. M., & Lo, Y. Y. 2009, ApJ, 691, 2
  * Kulsrud (2005) Kulsrud, R. M. 2005, Plasma Physics for Astrophysics (Princeton, NJ: Princeton Univ. Press)
  * Kuwabara & Ko. (2020) Kuwabara, T., & Ko, C. M. 2020, ApJ, 899, 72
  * Kuwabara et al. (2004) Kuwabara, T., Nakamura, K., & Ko, C. M. 2004, ApJ, 607, 828
  * Kuwabara & Ko (2006) Kuwabara, T., & Ko, C. M. 2006, ApJ, 636
  * Kuwabara & Ko (2015) Kuwabara, T., & Ko, C. M. 2015, ApJ, 798, 79
  * Li et al. (2015) Li, M., Ostriker, J. P., Cen, R., Bryan, G. L., & Naab, T. 2015, ApJ, 814, 4
  * Lo et al. (2011) Lo, Y. Y, Ko, C. M., & Wang, C. Y. 2011, Computer Physics Communications, 182, 177
  * Mao & Ostriker (2018) Mao, S. A., & Ostriker, E. C. 2018, ApJ, 854, 89
  * Parker (1958) Parker, E. N. 1958, ApJ, 128, 664
  * Parker (1963) Parker, E. N. 1963, Interplanetary Dynamical Processes (New York: Interscience Publ.)
  * Parker (1966) Parker, E. N. 1966, ApJ, 145, 811
  * Parker (1969) Parker, E. N. 1969, Space Sci. Rev., 9, 3
  * Ramzan et al. (2020) Ramzan, B., Ko, C.M.& Chernyshov, D.O. 2020, ApJ, 117, 905
  * Recchia (2020) Recchia, S. 2020 Int J Mod Phys D 29, 07
  * Recchia et al. (2016) Recchia, S., Blasi, P., & Morlino G., 2016, MNRAS, 462, 4227
  * Ruszkowski et al. (2017) Ruszkowski, M., Yang, H.-Y. K., & Zweibel, E. 2017, ApJ, 834, 208
  * Ryu et al. (2003) Ryu, D., Jongsoo, K., Seung S. H., & Jones, T. W. 2003, ApJ, 589, 1
  * Skilling (1975) Skilling, J. 1975, MNRAS, 172, 557
  * Wiener et al. (2017) Wiener, J., Pfrommer, C., & Oh, S. P. 2017, MNRAS, 467, 646
  * Yang et al. (2012) Yang, H.-Y. K., Ruszkowski, M., Ricker, P. M., Zweibel, E., & Lee, D. 2012, ApJ, 761, 2
  * Yu et al. (2020) Yu, B. P. B., Owen, E. R., Wu, K., & Ferreras, I. 2020, MNRAS, 492, 3179
  * Zhang (2018) Zhang, D. 2018, Galaxies, 6, 114
  * Zirakashvili et al. (1996) Zirakashvili, V. N., Breitschwerdt, D., Ptuskin, V. S., & Völk, H. J. 1996, A&A, 311, 113
  * Zweibel (2017) Zweibel, E. G. 2017, PhPl, 24, 055402

##  Appendix A Critical behaviors of three-fluid systems

In this appendix, we analyse the critical behaviour (or fixed points) of two
three-fluid systems. This is important for the study of transonic solutions
(i.e., solutions passing through fixed points or critical points).

###  A.1 Three-fluid system with cooling but without cosmic ray diffusion and
wave damping

We considered a three-fluid system (one-wave system with forward-propagating
wave) without cosmic ray diffusion and wave damping. The cosmic ray pressure
and wave pressure can be expressed in terms of flow speed by solving Equations
(3) and (3) (see Ramzan et al. 2020):

$P_{c}=A_{c}\left|\frac{\rho}{(1+\tilde{\psi}\sqrt{\rho})}\right|^{\gamma_{c}}\,,$
(23)
$P^{+}_{w}=\left[A^{+}_{w}-\frac{\gamma_{c}P_{c}}{2(\gamma_{c}-1)}\frac{\tilde{\psi}\sqrt{\rho})}{\left(1+\tilde{\psi}\sqrt{\rho})\right)}\right]\frac{\rho^{3/2}}{\left(1+\tilde{\psi}\sqrt{\rho})\right)^{2}}\,,$
(24)

where $\tilde{\psi}=\psi^{\prime}_{B}/(\sqrt{\mu_{0}}\psi^{\prime}_{m})$,
$\rho=\psi^{\prime}_{m}/U$. The governing set of equations can be written as a
three-dimensional autonomous system,

$\frac{d}{ds}\begin{bmatrix}\xi\\\ U\\\ \
P_{g}\end{bmatrix}=\begin{bmatrix}\mathcal{P}\\\ \mathcal{Q}\\\
\mathcal{R}\end{bmatrix},$ (25)

where

$\displaystyle\mathcal{P}=\left(1-M_{\rm eff}^{-2}\right),$ (26)
$\displaystyle\mathcal{Q}=\frac{1}{U}\left[g_{e}+\frac{(\gamma_{g}-1)\Gamma}{\psi^{\prime}_{m}}\right],$
(27)
$\displaystyle\mathcal{R}=-\frac{\gamma_{g}p_{g}}{U}\mathcal{Q}-\frac{(\gamma_{g}-1)\Gamma}{U}\mathcal{P}\,,$
(28)

where $M_{\rm eff}$= $U/a_{\rm eff}$, $a_{\rm eff}^{2}$ is given by Equation
(18). The fixed point or critical point of the system is determined by setting
the right hand side of the Equation (25) equal to zero. This can be achieved
by:

$\displaystyle\mathcal{P}=0,\quad\quad\mathcal{Q}=0,$ (29)

which gives a critical line in the three-dimensional ($\xi$,$U$,$P_{g}$)
space,

$\displaystyle U_{s}=\mu,$ (30)
$P_{gs}=\frac{1}{\gamma_{g}}\left[\psi^{\prime}_{m}U_{s}-\gamma_{c}P_{cs}\frac{(1+\frac{1}{2}\tilde{\psi}\sqrt{\rho_{s}})^{2}}{(1+\tilde{\psi}\sqrt{\rho_{s}})^{2}}-\frac{3}{2}P^{+}_{ws}\frac{(1+\frac{1}{3}\tilde{\psi}\sqrt{\rho_{s}})}{(1+\tilde{\psi}\sqrt{\rho_{s}})}\right]\,,$
(31)
$\xi_{s}=\sqrt{\frac{GM\psi^{\prime}_{m}}{(\gamma_{g}-1)\Gamma_{s}}}-a_{K},$
(32)

where $\mu$ is the parametric variable of the critical line and
$\Gamma_{s}=\Gamma_{s}(\rho_{s},P_{gs})$. For illustration, we take
$\Gamma=\delta_{c}\rho^{2}(P_{g}/\rho)^{1/2}$ and show the result in Figure 6.
The critical line is the long-dashed pink line. Linear analysis show that the
points on the critical line are saddle-points except for those close to
$P_{g}=0$. In the figure, we show transonic solutions from three points on the
critical line. For each point, there are two transonic solutions: one from
subsonic to supersonic (solid line) and the other from supersonic to subsonic
(dashed line).

Figure 6:  Critical line and three sets of transonic solutions of a three-
fluid system without cosmic ray diffusion. In this figure, $\gamma_{g}=5/3$,
$\gamma_{c}=4/3$, $GM=1$, $a=1$, $\psi^{\prime}_{B}=1/2$,
$\psi^{\prime}_{m}=1$, $\delta_{c}=0.1$, $A_{c}=0.5$, $A_{w}^{+}=0.25$.
Magenta long-dashed line is the critical line. Solid lines: subsonic-
supersonic transitions. Dashed lines: supersonic-subsonic transitions.

###  A.2 Three-fluid system with cosmic ray diffusion but without cooling and
wave damping

We also considered a three-fluid system (with forward-propagating wave) in
which cosmic ray diffusion is included but cooling and wave damping are
neglected. The diffusion coefficient is inversely proportional to the wave
pressure. In addition to magnetic flux and mass flux, the total energy flux
and the wave-action integral are constant in this system:

$\displaystyle F^{\prime}_{\rm tot}=\frac{1}{2}U^{2}+F^{\prime}_{\rm
g}+F^{\prime}_{\rm c}+F^{\prime+}_{w}+\Psi_{e}\,,$ (33)
$\displaystyle\mathcal{W}_{A}^{\prime}=\left[F^{\prime}_{\rm
c}+\frac{2P^{+}_{w}U}{\psi^{\prime}_{m}}\frac{(1+\tilde{\psi}\sqrt{\rho})^{2}}{\tilde{\psi}\sqrt{\rho}}\right]\,,$
(34)

$F^{\prime}_{g}=F_{g}/\psi^{\prime}_{m}$,
$F^{\prime}_{c}=F_{c}/\psi^{\prime}_{m}$,
$F^{\prime+}_{w}=F_{w}^{+}/\psi^{\prime}_{m}$ and $F_{g}$, $F_{c}$ and
$F_{w}^{+}$ are given by Equations (7), (2), and (9). Here $\Psi_{e}$ is the
external gravitational potential and:

$\rho U=\psi^{\prime}_{m},\quad
P_{g}=A_{g}\rho^{\gamma_{g}},\quad\kappa=\eta/3P^{+}_{w}.$ (35)

Solving Equations (33) and (34) for $P^{+}_{w}$ and $F_{c}^{\prime}$,

$\displaystyle P^{+}_{w}=$
$\displaystyle\,\frac{\tilde{\psi}\rho^{3/2}}{2(1+\frac{1}{2}\tilde{\psi}\sqrt{\rho})}\left[\frac{1}{2}U^{2}+\frac{\gamma_{g}P_{g}}{(\gamma_{g}-1)\rho}+\Psi_{e}-F^{\prime}_{\rm
tot}+\mathcal{W}_{A}^{\prime}\right]\,,$ (36) $\displaystyle F_{c}^{\prime}=$
$\displaystyle\frac{(1+\tilde{\psi}\sqrt{\rho})^{2}}{(1+\frac{1}{2}\tilde{\psi}\sqrt{\rho})}\left[F^{\prime}_{\rm
tot}-\frac{1}{2}U^{2}-\frac{\gamma_{g}P_{g}}{(\gamma_{g}-1)\rho}-\Psi_{e}\right]$
$\displaystyle\quad\quad-\frac{3\tilde{\psi}\sqrt{\rho}(1+\frac{2}{3}\tilde{\psi}\sqrt{\rho})}{2(1+\frac{1}{2}\tilde{\psi}\sqrt{\rho})}\mathcal{W}_{A}^{\prime}\
\,,$ (37)

The governing set of equations can be written as a set of three-dimensional
autonomous systems:

$\frac{d}{ds}\begin{bmatrix}\xi\\\ U\\\ \
P_{c}\end{bmatrix}=\begin{bmatrix}\mathcal{P}\\\ \mathcal{Q}\\\
\mathcal{R}\end{bmatrix},$ (38)

where

$\displaystyle\mathcal{P}=\left(1-M_{\rm eff}^{-2}\right),$ (39)
$\displaystyle\mathcal{Q}=$ $\displaystyle\,\frac{g_{e}}{U}-\mathcal{D}\,,$
(40) $\displaystyle\mathcal{R}=$
$\displaystyle\,\frac{(1+\tilde{\psi}\sqrt{\rho})}{(1+\frac{1}{2}\tilde{\psi}\sqrt{\rho})}\psi^{\prime}_{m}\mathcal{DP}\,,$
(41) $\displaystyle\mathcal{D}=$
$\displaystyle\frac{(\gamma_{c}-1)}{\kappa}(1+\tilde{\psi}\sqrt{\rho})$
$\displaystyle*\left[\frac{1}{2}\right.U^{2}+\frac{\gamma_{g}P_{g}}{(\gamma_{g}-1)\rho}+\frac{\gamma_{c}P_{c}}{(\gamma_{c}-1)\rho}\frac{(1+\frac{1}{2}\tilde{\psi}\sqrt{\rho})}{(1+\tilde{\psi}\sqrt{\rho})}$
$\displaystyle+\Psi_{e}-F^{\prime}_{\rm
tot}\left.+\frac{3\tilde{\psi}\sqrt{\rho}(1+\frac{2}{3}\tilde{\psi}\sqrt{\rho})}{2(1+\tilde{\psi}\sqrt{\rho})^{2}}\mathcal{W}_{A}^{\prime}\right]\,,$
(42)

where $M_{\rm eff}$= $U/a_{\rm eff}$, $a_{\rm eff}^{2}$ is given by Equation
(12) with $P_{w}^{-}=0$. The fixed point or critical point of the system is
determined by setting the right hand side of Equation (38) equal to zero. This
can be achieved via:

$\displaystyle\mathcal{P}=0,\quad\quad\mathcal{Q}=0,$ (43)

which gives a critical line in three dimensional ($\xi$,U,$P_{c}$) space,

$\displaystyle\Psi_{es}=$
$\displaystyle\left[F^{\prime}_{\rm{tot}}-\mathcal{W}_{A}^{\prime}-\frac{1}{2}U^{2}_{s}-\frac{\gamma_{g}P_{gs}}{(\gamma_{g}-1)\rho_{s}}\right]$
$\displaystyle+\frac{4(1+\frac{1}{2}\tilde{\psi}\sqrt{\rho_{s}})(1+\tilde{\psi}\sqrt{\rho_{s}})}{3\tilde{\psi}\sqrt{\rho_{s}}(1+\frac{1}{3}\tilde{\psi}\sqrt{\rho_{s}})}\left(U^{2}_{s}-\frac{\gamma_{g}P_{gs}}{\rho_{s}}\right),$
(44) $\displaystyle U_{s}=\mu,$ (45)
$\displaystyle\frac{\gamma_{c}P_{cs}}{(\gamma_{c}-1)\rho_{s}}=$
$\displaystyle\frac{\kappa_{s}g_{es}}{(\gamma_{c}-1)U_{s}(1+\frac{1}{2}\tilde{\psi}\sqrt{\rho_{s}})}$
$\displaystyle-\frac{4(1+\tilde{\psi}\sqrt{\rho_{s}})^{2}}{3\tilde{\psi}\sqrt{\rho_{s}}(1+\frac{1}{3}\tilde{\psi}\sqrt{\rho_{s}})}\left(U^{2}_{s}-\frac{\gamma_{g}P_{gs}}{\rho_{s}}\right)$
$\displaystyle+\frac{\mathcal{W}_{A}^{\prime}}{(1+\tilde{\psi}\sqrt{\rho_{s}})}.$
(46)

where $\Psi_{es}=\Psi_{e}(\xi_{es})$, $g_{es}=g_{e}(\xi_{es})$,
$\kappa_{s}=\kappa(P^{+}_{ws})$ and $\rho_{s}$, $P_{gs}$ and $P^{+}_{ws}$ are
given by Equations (35) and (36); also, $\mu$ is the parametric variable of
the critical line.

##  Appendix B Summary of three-fluid outflows

In this appendix, we discuss some general properties of the three-fluid system
with a forward-propagating wave and cosmic-ray diffusion. In the realm of
physically allowable solutions starting in the subsonic region at the base of
potential, the profiles of $U$, $P_{c}$, and $P_{g}$ can be categorized into
two regimes (Ramzan et al. 2020): regime A (called cosmic ray accompanied
outflows) and regime B (quasi-thermal outflows). In regime A, the value of $U$
increases ($P_{c}$ and $P_{g}$ decreasing) gradually to a constant value at
large distances (e.g., Figure 1). Cosmic rays are coupled to the thermal gas
at all distances (i.e., $P_{w}^{+}$ remains finite all the way). In regime B,
values of $U$ and $P_{c}$ have a general trend of increasing towards a finite
value at large distances ($P_{g}$ decreases). The profile of $U$ and $P_{c}$
has a characteristic rapid increase or a “jump” at some distance from the base
($P_{g}$ has a rapid decrease or a “drop,” e.g., Figure 2). Also, $P_{w}^{+}$
decreases monotonically and has a “drop” at the same location and approaches
zero afterward. Cosmic rays are decoupled from the thermal gas after the
“jump” and $P_{c}$ becomes constant afterward. Wave damping and cooling are
not necessary conditions for the existence of the two regimes, although they
may affect parameters defining the regimes. In regime A, the set of boundary
conditions at the base that gives a physically allowable solution is
“isolated” in the sense that any slight change in one (and only one) parameter
will push the system into the unphysical domain (however, if more than one
parameter can be changed, then it is possible to get into some other
physically allowable domain). In contrast, in regime B, a small change in one
(and only one) parameter of a physically allowable set may also give physical
solutions. This can be understood as follows.

To simplify the argument, let us consider a basic case (system without wave
damping and cooling). In this case, only the first two terms on the right hand
side of Equation (11) remain (with $e_{+}=1$ and $e_{-}=0$), and the second
and third terms on the right hand side of Equation (16) vanish. Equation (16)
can be viewed as an equation for the cosmic ray diffusive flux $D_{c}=\kappa
dE_{c}/d\xi$ ($\kappa=\eta/3/P_{w}^{+}$). When $M_{\rm eff}<1$ (“subsonic”
region), Equation (11) and Equation (16) form a positive feedback loop in the
following sense. If both $dU/d\xi$ and $D_{c}$ are negative (positive), then
both will be more negative (positive) downstream. There are two possible
routes, route (a) and route (b), for avoiding divergence and ending up with a
physically allowable solution.

Let the solution start at the base with $D_{c}<0$ and $dU/d\xi>0$. $D_{c}$
then increases while $dU/d\xi$ decreases. For one particular set of boundary
values, $D_{c}$ becomes less and less negative and approaches to zero at large
distances (and $dU/d\xi$ approaches zero from the positive side). This is
route (a).

If a solution starts with the same boundary values as route (a) except that
$D_{c}$ is (slightly) smaller, then $dU/d\xi$ will turn negative before
$D_{c}$ goes to zero. The positive feedback loop activates and the solution
will go to the unphysical domain ($U$ and $P_{c}$ become negative). However,
if a solution starts with the same boundary values as route (a) except that
$D_{c}$ is larger, the situation will be different.

Supposing the solution starts from $D_{c}<0$ and $dU/d\xi>0$, $D_{c}$
increases and becomes positive at some finite $\xi$. The positive feedback
loop kicks in and both $D_{c}$ and $dU/d\xi$ increase rapidly. A positive
$D_{c}$, namely, a positive cosmic-ray pressure gradient, will cause
$P_{w}^{+}$ to diminish via the cosmic-ray streaming instability. Cosmic rays
will be decoupled from the thermal gas and the outflow behaves like a thermal
outflow. If the velocity and thermal pressure at that location belong to the
physically allowable region of a pure thermal outflow, then a physically
allowable solution is obtained: this is route (b). The set of parameters that
allows for route (b) belonging to regime B (quasi-thermal). If the velocity
and thermal pressure after the “rapid change” do not belonging to any
physically allowable region of a pure thermal outflow, then the set of
parameters only allows for route (a) and it belongs to regime A mentioned
above.

Once again, the characteristics of a quasi-thermal outflow is a “jump” in
velocity and cosmic-ray pressure (and a “drop” in thermal pressure) at the
location where $P_{w}^{+}$ diminishes towards zero. A small change in $D_{c}$
at the base will change the location of the “jump.” Smaller $D_{c}$ (more
negative) pushes the “jump” to larger distances. We point out that there is a
minimum $D_{c}$ such that the “jump” occurs at “infinity” (e.g., the left
column of Figure 2). This particular solution coincides with route (a). The
solutions become unphysical for value of $D_{c}$ that are smaller than the
minimum.