# Rapid CO gas dispersal from NO Lup’s class III circumstellar disc

J. B. Lovell1, G. M. Kennedy2,3, S. Marino1, M. C. Wyatt1, M. Ansdell4, M.
Kama5,1, C. F. Manara6, L. Matrà7, G. Rosotti8, M. Tazzari1, L. Testi6,9, J.
P. Williams10
1Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge,
CB3 0HA, UK
2Department of Physics, University of Warwick, Coventry, CV4 7AL, UK
3Centre for Exoplanets and Habitability, University of Warwick, Gibbet Hill
Road, Coventry CV4 7AL, UK
4National Aeronautics and Space Administration Headquarters, 300 E Street SW,
Washington DC 20546, USA
5Tartu Observatory, University of Tartu, 61602 Tőravere, Estonia
6European Southern Observatory, Karl-Schwarzschild-Strasse 2, 85748, Garching
bei München, Germany
7School of Physics, National University of Ireland Galway, University Road,
Galway, Ireland
8Leiden Observatory, Leiden University, P.O. Box 9513, NL-2300 RA Leiden, the
Netherlands
9INAF - Osservatorio Astrofisico di Arcetri, L.go E. Fermi 5, 50125, Firenze,
Italy
10Institute for Astronomy, University of Hawai’i at Mänoa, Honolulu, HI, USA
E-mail<EMAIL_ADDRESS>
(Accepted 2020 November 17. Received 2020 November 17; in original form 2020
October 16)

###### Abstract

We observed the K7 class III star NO Lup in an ALMA survey of the 1-3 Myr
Lupus association and detected circumstellar dust and CO gas. Here we show
that the J = 3-2 CO emission is both spectrally and spatially resolved, with a
broad velocity width ${\sim}19$ km s-1 for its resolved size ${\sim}1\arcsec$
(${\sim}130$ au). We model the gas emission as a Keplerian disc, finding
consistency, but only with a central mass of ${\sim}11M_{\odot}$, which is
implausible given its spectral type and X-Shooter spectrum. A good fit to the
data can also be found by modelling the CO emission as outflowing gas with a
radial velocity ${\sim}22$ km s-1. We interpret NO Lup’s CO emission as the
first imaged class III circumstellar disc with outflowing gas. We conclude
that the CO is continually replenished, but cannot say if this is from the
break-up of icy planetesimals or from the last remnants of the protoplanetary
disc. We suggest further work to explore the origin of this CO, and its higher
than expected velocity in comparison to photoevaporative models.

###### keywords:

circumstellar matter - planetary systems - submillimetre: planetary systems.

††pubyear: 2020††pagerange: Rapid CO gas dispersal from NO Lup’s class III
circumstellar disc–A

##  1 Introduction

Stars are born with protoplanetary discs containing large quantities of
primordial gas and dust which persist for several Myr before dispersing on
rapid ${\sim}$0.1 Myr timescales (Ercolano & Pascucci, 2017). Circumstellar
discs are also seen around older stars (${\gtrapprox}$10 Myr), known as debris
discs, where the dust and gas is inferred to be secondary, created in the
destruction of planetesimals that must be volatile-rich to replenish the gas
(Wyatt, 2008; Dent et al., 2014; Marino et al., 2016; Moór et al., 2017; Matrà
et al., 2017; Kral et al., 2017). The transition between the two types of disc
is not well understood, but is thought to involve a combination of gas being
accreted onto the star, or being expelled from the system by disc winds driven
by photoevaporation or magneto-hydrodynamics (MHD), as well as planet
formation processes (Williams & Cieza, 2011; Wyatt et al., 2015; Lesur, 2020).
Class III stars are those in star forming regions for which infrared emission
shows an absence of hot dust suggesting that the star’s protoplanetary disc
has either recently dispersed or is in the process of dispersal (Adams et al.,
1987). However, these stars usually have limited constraints on the presence
of cold dust or gas with which to constrain their nature.

A recent ALMA survey of class III stars in Lupus found dust in several systems
with dust masses orders of magnitude lower than protoplanetary disc levels and
consistent with originating in the progenitors of debris discs seen at later
ages, though this does not rule out the possibility that this dust is a
remnant of the protoplanetary disc (Lovell et al., 2020). For one of these
class III stars, NO Lup, J = 3-2 CO gas was also detected, allowing a more
thorough assessment of the nature of its circumstellar environment. Comparison
of the inferred mass of CO with that of the 10 Myr old M star TWA 7, for which
the level was consistent with secondary production (i.e.,
$0.8{-}80{\times}10^{-6}M_{\oplus}$, see Matrà et al., 2019), shows that the
CO could indeed be secondary, though as for the dust interpretation, such
plausibility does not preclude the possibility that the CO is primordial.

This paper presents a detailed analysis of the CO detected towards NO Lup,
using spectral and spatial data not reported in Lovell et al. (2020),
constraining the kinematic structure of the gas disc. While circumstellar gas
is usually seen in a Keplerian disc, other gas morphologies are possible, such
as an outflowing component which could be a feature of either a depleting
primordial disc (Haworth & Owen, 2020) or a debris disc. In $\S$2 we provide a
summary of the NO Lup system, and introduce new X-Shooter observations. In
$\S$3 we extend the analysis of the CO gas by showing the problem with
modelling this as a Keplerian disc, and show better consistency in $\S$4 by
modelling this with an outward radial velocity. We interpret our results in
$\S$5, and conclude in $\S$6.

##  2 NO Lup in Context

NO Lup (2MASS J160311.8-323920) is located in Cloud I of Lupus at
${\alpha}{=}$16:03:11.812, ${\delta}{=}$-32:39:20.31 (J2000) at a distance of
$133.7{\pm}0.7$ pc (Gaia Collaboration et al., 2018). NO Lup has previously
been classified spectrally as non-accreting (Cieza et al., 2013), and as a K7
star (Krautter et al., 1997). The latter is consistent with its Gaia DR2
temperature, $T_{\rm{eff}}{=}$3994$K$, and stellar luminosity,
$L_{\star}{=}0.287L_{\odot}$ and its well constrained spectral energy
distribution (given observations with WISE, Spitzer, and 2MASS, which found
the blackbody planetesimal belt radius as $R_{\rm{BB}}{=}3.2{\pm}0.3$ au, see
Lovell et al., 2020). Analysing the X-Shooter spectrum of NO Lup (project
093.C-0506 A), we confirm both spectral and accretion analyses. By reducing
this data using the standard pipeline version 3.5.0 (Modigliani et al., 2010),
we show in Fig. 1 (top and middle) the full and zoomed-in X-Shooter spectra
for NO Lup, which are consistent with the spectral features of the well
characterised class III K7 star, SO879. In the lower panel we show the
H${\alpha}$ emission line with an
EW=$-2.80{\pm}0.15{$\mathrm{\SIUnitSymbolAngstrom}$}$, which we find to be
centred on a radial velocity (RV) of $-3.5{\pm}2.0$ km s-1, consistent with
the Gaia DR2 RV of $-1.93{\pm}4.08$ km s-1, and the average RV of Lupus stars,
$\rm{RV_{Lup}}=2.8{\pm}4.2\rm{km\,s^{-1}}$ (Frasca et al., 2017). This EW is
consistent with non-accretion and only slightly higher than, but also
consistent with, the line width of non-accreting SO879. Hardy et al. (2015)
estimated the stellar mass of NO Lup to be $M_{\star}{=}0.7M_{\odot}$. With
the models of Siess et al. (2000) and Baraffe et al. (2015), and the Gaia DR2
$L_{\star}$ and Teff, we estimated $M_{\star}$ between 0.7-0.8$M_{\odot}$,
consistent with the literature. Although NO Lup has significant emission above
the photosphere at 12 and 24$\mu$m, these excesses are small, resulting in
mid-IR spectral slopes steeper than those of protoplanetary discs. Unresolved
continuum emission was detected by ALMA, implying a dust mass of
$0.036{\pm}0.007M_{\oplus}$ and disc radius <56 au (Lovell et al., 2020).
These observations also detect CO J=3-2 line emission, with
$F_{\rm{CO}}{=}0.29{\pm}0.07$ Jy km s-1, and a width of ${\sim}19$ km s-1
(between -11.0 to +8.1 $\rm{km}\,\rm{s}^{-1}$), consistent with being centred
on the stellar RV discussed above. Assuming this emission is optically thin
and in LTE at $T{=}50K$, this flux corresponds to a CO gas mass of
$4.9{\pm}1.1{\times}10^{-5}M_{\oplus}$, which for an ISM CO abundance implies
a gas-to-dust ratio of $1.0{\pm}0.4$ (Lovell et al., 2020). We note that a
wider temperature range of $20{-}100$ K, consistent with the range of
literature gas temperatures, could change this gas mass by at most a factor of
${\sim}$2 (see equations 2 and 8 of Matrà et al., 2017).

Figure 1: X-Shooter spectra for NO Lup and SO879. Top: full X-Shooter range.
Middle: zoomed-in 690-720 nm range. Bottom: H${\alpha}$ line. SO879’s flux is
scaled to NO Lup’s median source flux between 698-702 nm, and is wavelength
corrected to have a common RV with NO Lup based on SO879’s DR2 RV of $30.10\pm
0.21$ km s-1 (Gaia Collaboration et al., 2018). Figure 2: Averaged channel map
(4 ${\sim}$0.4 km s-1 channels per plot), showing the $\pm 3$ and $5\sigma$
emission significance. In all plots North is up, East is left, contour lines
represent $\pm 3$ and $\pm 5\sigma$, 100 au scale bars are shown in the lower-
left of the bottom panel plots, the synthesised beam is shown in the bottom-
right plot, and the major axis ticks are $1^{\prime\prime}$ wide.

##  3 The Problem with a Keplerian Disc

Assuming the gas is in Keplerian rotation, the width of the CO line can be
used to estimate the spatial extent of the emission. For example, if the star
has a mass ${\sim}0.7M_{\odot}$ (see $\S$2), the line width indicates a radius
of ${\sim}10$ au (or smaller if the disc is not edge on), well below the
spatial resolution limit of our measurements. However, we find that the CO
emission is spatially resolved. Fig. 2 shows that the centroid of the CO
emission transitions from the South-West to the North-East when binned in
channels with decreasing RV, as expected for emission that originates in a
disc inclined to the plane of the sky. To resolve this emission given a beam
of ($0.94\times 0.82$)$\arcsec$ (at position angle $70^{\circ}$), the disc
emission would have to have a radius ${\gtrapprox}\,50$ $\rm{au}$. This
however would then be inconsistent with the previously inferred ${\sim}10$ au
radius for the assumed stellar mass. This conclusion can also be illustrated
with the position-velocity (PV) diagram in Fig. 3 (see left-hand plot 1),
produced with a 36∘ position angle and 2$\arcsec$ slit width, which shows how
the offset varies with radial velocity. This agrees with the previous analysis
that the velocities are unexpectedly high at a large separation (nearly
$10\,\rm{km}\,\rm{s}^{-1}$ at ${\sim}\,0.5\arcsec{=}66$ au). The curves and
radial extension line on this same plot show that for gas to be in a Keplerian
orbit the stellar mass would have to be ${\sim}11M_{\odot}$, i.e., much higher
than $0.7M_{\odot}$.

We use a simple model to explore this Keplerian disc interpretation by
assuming the gas is in a Keplerian orbit with the stellar mass of NO Lup left
as a free parameter, $M_{\star}$. This Disc Model computes the density of
emitting CO for each pixel in a cube that has RA, Dec, and line of sight, $z$,
as axes. The orientation, position and RV of the disc in the cube are set by
the ascending node, $\Omega$ (i.e., the PA of the disc major axis), the
inclination, $i$, a phase centre offset in RA and Dec (xoff and yoff), the
systemic velocity, $v_{\rm{sys}}$, and the fixed distance to NO Lup of 133.7
pc (see $\S$2). The model assumes that the CO density can be modelled with a
total flux, $F$, and a power-law radial profile, defined between $r_{\rm{in}}$
and $r_{\rm{out}}$ with index $p$, for which the volume density goes as
$r^{p-1}$. The model assumes a Gaussian scale height with a fixed aspect
ratio, $\sigma_{h}=0.05$, meaning that the surface density goes as
$\Sigma{\propto}r^{p}$. We fitted the visibilities of channels -15.1 to
${+}6.9$ km s-1 from the CO line (as done in Kennedy et al., 2019), using the
emcee package (Foreman-Mackey et al., 2013), with 40 walkers and 2000 steps to
achieve convergence. The fit gives reasonable results, showing no residuals in
the PV diagram ${>3}\sigma$ (see Fig. 3, plots 2 and 3), with the best-fit
model values shown in Table 1. We note that the best-fit
$v_{\rm{sys}}{=}{-3.8}$ km s-1 is consistent with that in $\S$2, and that this
model did not constrain $p$, which uniformly varied between the imposed limits
of -2 to 1. A similar result was found with a Gaussian radial distribution.
This fitting procedure therefore comes to the same conclusion as shown earlier
in requiring a stellar mass $M_{\star}{=}11.1{\pm}1.9M_{\odot}$ that is
significantly greater than expected. Since the stellar luminosity rules out
${\sim}10$ K-type stars within 10s of au of NO Lup, and the X-Shooter spectra
rules out NO Lup being a mis-classified star, we conclude that we cannot
consistently interpret the CO emission as originating from a Keplerian disc.

Table 1: Best-fit model parameters. D: Disc Model, and O: Outflow Model. Mod | $\Omega$ | $i$ | $r_{\rm{in}}$ | $r_{\rm{out}}$ | $M_{\star}$ | $v_{r}$
---|---|---|---|---|---|---
| deg | deg | au | au | $M_{\odot}$ | km s-1
D | $36{\pm}3$ | $42{\pm}5$ | $35{\pm}8$ | $90{\pm}5$ | $11.1{\pm}1.9$ | -
O | $115{\pm}3$ | $20{\pm}4$ | $20{\pm}4$ | $120{\pm}14$ | $0.70{\pm}0.05$ | $22{\pm}4$

##  4 An Outflowing Gas Interpretation

Having ruled out the Keplerian disc interpretation, we explore the possibility
that the gas emission is dominated by a radially outflowing velocity
component. We define an Outflow Model with a stellar mass fixed at
$0.7M_{\odot}$ (which sets the azimuthal velocity), where the CO emission
extends between $r_{\rm{in}}$ and $r_{\rm{out}}$ with a radial velocity,
$v_{r}$, which models the CO gas outflow. This is a variant of the Disc Model
of $\S$3, with the same disc parameters and geometry, except that the gas also
has an additional radial velocity component and the vertical density
distribution is modelled as uniform, with fixed lower and upper edges, such
that a vertical cross-section through the disc looks like a wedge with an
opening angle $\delta h$. Allowing for a larger scale height is a rough
approximation to disc wind models, where material flows both vertically off
the disc and radially outwards. We fitted the visibilities as described above
for the Disc Model. We again find reasonable results (see Fig. 3, plots 4 and
5), and quote the best-fit model parameters in Table 1, also finding
$v_{\rm{sys}}{=}-3.8\rm{km}\,\rm{s}^{-1}$, and an unconstrained $p$, which we
therefore fixed as $p{=}-1$111$p{=}-1$ is the value expected for mass
conservation if the CO is being created at the inner disc edge.. The Outflow
Model finds $v_{r}{\approx}22$ km s-1, which is higher than the measured
velocity from the CO line width and thus dominates the azimuthal velocity for
the modelled $r_{\rm{in}}$ and $r_{\rm{out}}$. The Outflow Model finds a large
opening angle of ${\delta}h{=}0.3{\pm}0.1$; the disc inclination is required
to be low to reproduce the spatial extent, thus the scale height is large to
maximise the radial component of the velocity, and $v_{r}$ is significantly
larger than the minimum possible value of ${\sim}19/2$ km s-1. Though the flow
in the model is mostly radial, the large scale height can be thought of as
approximating the significant vertical component present in outflow models.
Despite the simplicity of assuming the gas to have an additional radial
velocity component, this model shows that the observations are broadly
consistent with a scenario in which the gas is radially outflowing.

Figure 3: Plot 1: PV diagram showing the distribution of velocities present at
different offsets along a slit width of 2$\arcsec$ and position angle of 36∘.
The contours show 3 and $5\sigma$ emission. The curves demonstrate the maximum
allowed radial velocities for Keplerian motion around a $11M_{\odot}$ star
(green solid) and a $0.7M_{\odot}$ star (blue dashed), and the radial velocity
expected for a disc with an 80 au radius (black dot-dash slope). The lower-
left thumbnail shows the Moment-0 map (box extent ${\sim}300$ au, with 3, 6,
and 9${\sigma}$ contours). Plots 2 and 4: Model PV diagrams using wider slits
covering the full model extent, however we note that the emission beyond
$2\arcsec$ is negligible. The thumbnails show the model Moment-0 maps
(${\sim}300$ au box widths), and PAs (see Table 1). Plots 3 and 5: Residual PV
diagrams, with $\pm 2$ and $\pm 3\sigma$ residuals shown in blue/orange, and
the same contours as the data (black dashed).

##  5 First Imaged Class III Gas Dispersal

Whilst we model the CO emission over all the channels in which it is detected,
for simplicity we have plotted the two models of $\S$3 and 4 in Fig. 3 as PV
diagrams, with their corresponding moment-0 maps as thumbnails. These show
that whereas the models have similar PV diagrams, the position angles of their
disc midplanes are ${\sim}90^{\circ}$ different, discussed further in Appendix
A. At low spatial resolution, Keplerian disc models and high radial velocity
outflowing models (i.e., those in which the radial velocity dominates the
azimuthal velocity) can have indistinguishable PV diagrams. This degeneracy
can be broken with higher resolution imaging to measure the PA of the
continuum emission which would confirm the radial (or azimuthal) nature of
observed CO velocities.

The simple Outflow Model presented in $\S$4 showed that the observations can
be fitted with a constant ${\sim}22$ km s-1 radial velocity component. While
it may be possible to explain the observations with a lower outflow velocity,
for example with a more detailed model of the gas kinematics, any gas flow
must be ${\gtrapprox}10$ km s-1, given the spectrum line width of ${\sim}19$
km s-1. This is high compared to the photoevaporative models of Haworth & Owen
(2020), in which the outflowing wind velocity is ${\sim}3$ km s-1. In fact,
this high velocity rules out a scenario in which the gas is a pure
photoevaporative wind, since such a wind velocity should be set by the sound
speed, and so to reach even ${\sim}10$ km s-1 would require the CO gas to be
at $T{\sim}10,000$ K, which exceeds the CO thermal dissociation temperature.
To accelerate CO gas to such high speeds, additional forces are required
(e.g., such as MHD driven winds, see Lesur, 2020, and references therein).
Winds thought to be driven by the disc magnetic field have been observed
towards less evolved class II stars with similarly high velocities (see
Pontoppidan et al., 2011). Theoretically, magnetically driven disc winds can
be launched with velocities ${\sim}$several factors of the Keplerian velocity
at 10s of au (Bai, 2016), which at the location of our disc inner edge
($r_{\rm{in}}{=}20{\pm}4$ au) is consistent with our modelled outflow
velocity. Thus the 10s of km s-1 velocities inferred here are plausible via
this mechanism. Though we are not aware of theoretical works that have looked
into the launching of a magnetised disc wind in the conditions of a debris
disc, we note that 4-12 km s-1 dust outflows were detected towards AU Mic by
Boccaletti et al. (2018), and it is worth exploring whether these are launched
by a similar mechanism. Interpretation of NO Lup’s CO emission may therefore
require modelling a 3D velocity field and MHD driving, in addition to CO
photodissociation and shielding (which we discuss below, Kral et al., 2019;
Marino et al., 2020).

While we have discussed the parameter $v_{r}$ in the Outflow Model as an
outflowing velocity, the fit to the observations is insensitive to the sign of
$v_{r}$, thus we next explore whether the gas is more likely to be infalling
or outflowing. To do so, we compare the observed CO mass loss rate with the
upper limit on the CO accretion rate. Assuming
$\dot{M}_{\rm{CO}}{=}M_{\rm{CO}}.v_{r}$/$R$, we find
$\dot{M}_{\rm{CO}}{\sim}3M_{\oplus}$ Myr-1, for
$M_{\rm{CO}}{\sim}5{\times}10^{-5}M_{\oplus}$, at $v_{r}{\sim}22$ km s-1, from
the mean disc radius $R{=}$($r_{\rm{in}}$+$r_{\rm{out}}$)/2${\sim}70$ au.
Since the estimate of $M_{\rm{CO}}$ assumes it is in LTE and that the gas is
optically thin, this mass loss rate is a lower limit. Next, by fitting a
polynomial to the continuum near the H${\alpha}$ line (with the EW stated in
$\S$2), we obtain a continuum-subtracted line flux of
${\sim}2.9{\times}10^{-15}$ erg s-1 cm-2, from which we estimate the accretion
luminosity as $\log{L_{\rm{acc}}/L_{\odot}}{\sim}{-3.04}$, using the empirical
relation between line luminosity and accretion luminosity of Alcalá et al.
(2017). With $L_{\rm{acc}}$, $M_{\star}$ and $L_{\star}$ from $\S$2, and a
stellar radius of $R_{\star}{\sim}1.3R_{\odot}$ (Gaia Collaboration et al.,
2018), we find a $3\sigma$ upper limit CO mass accretion rate of
<0.1$M_{\oplus}$ Myr-1, for an ISM H2/CO abundance ratio of 10,000 (i.e.,
consistent with the ratio in primordial gas). This is consistent with other
non-accreting class III stars (see Manara et al., 2013), however it is more
than 30 times lower than the inferred CO mass loss rate, indicating that if
primordial the CO gas cannot be inflowing. While this cannot rule out
inflowing gas for lower H2/CO ratios (e.g., if the gas is produced in a
secondary scenario), we rule out inflowing secondary gas later. However, the
upper limit on the CO mass accretion rate may provide an important constraint
on models for the outflow that require an inflowing component to conserve
angular momentum.

To explore the origin of the gas, we first note that CO at 10-20 km s-1 (i.e.,
2-4 au yr-1) would travel ${\gtrapprox}$200 au in 100 yr, i.e., over the CO
photodissociation time (Visser et al., 2009), although if well shielded the CO
could survive much longer. This may suggest that the gas was formed during the
recent break-up of a massive planetesimal, however such an event would be both
rare and leave a distinctive asymmetry in the gas distribution (Jackson et
al., 2014) which we do not observe. Rather, the measured ${\sim}130$ au extent
of the CO gas could still be consistent with a $>$100 yr lifetime, as at
larger distances from the star cooler CO is less collisionally excited, and so
difficult to detect (Matrà et al., 2015). Thus, comparing the extent, velocity
and gas lifetime does not lead to strong constraints, though this may not be
the case for neutral CI gas, which is likely a better probe of outflowing gas
at larger distances, as suggested by Haworth & Owen (2020). Given typical
stellar/disc timescales of 0.1-1 Myr, the high velocity however suggests that
the CO must be continuously replenished, as it is extremely unlikely that the
CO was imaged within 100 yr after a single or final production event.

The CO reservoir that replenishes the gas may either be in gaseous form (a
protoplanetary disc remnant) or in solid form (icy planetesimals in a debris
disc). Protoplanetary discs are expected to disperse on ${\sim}$100 kyr
timescales (e.g. Ercolano & Pascucci, 2017). While NO Lup’s SED suggests that
it has already lost its protoplanetary disc, the dust may have dispersed
first, leaving a primordial CO reservoir (e.g., as suggested by Owen &
Kollmeier, 2019). The detection of CO gas towards 1/30 of the 2 Myr old Lupus
class III stars in Lovell et al. (2020) implies plausible dispersal timescales
for such primordial gas remnants of ${\sim}$70 kyr. A potential problem with
this, however, is that no CO reservoir with a Keplerian disc signature is
present in our observations. For such a CO reservoir to go undetected its
surface brightness, which in the optically thick limit is
$I_{\nu}{=}B_{\nu}(T)A/A_{\rm beam}$ (where $A$ is the CO emission area, and
$A_{\rm beam}$ is the beam area), should be below the 3$\sigma$ noise level of
${\sim}9$ mJy beam-1 (see Lovell et al., 2020). If CO emission fills the beam
($A{=}A_{\rm beam}$), then this upper limit implies a temperature below
${\sim}$4 K, significantly below the CO sublimation temperature and thus
unlikely. For a more reasonable temperature of 50 K, the CO emitting area
would need to be ${\sim}$150 times smaller than the beam, and thus a ring at
$r_{\rm{in}}{\sim}20$ au must have a width narrower than 0.01 au, which is
also unlikely. This argues against an optically thick ring of CO being the gas
source, though further high-resolution imaging is required to definitively
conclude this. A more plausible explanation may be that the wind is
replenished by a reservoir of CO in icy planetesimals in the
$R_{\rm{BB}}{\sim}3$ au belt. Since blackbody disc radii estimates can be
${\gtrapprox}5{\times}$ smaller than physical planetesimal belt radii (see
Eq.8 of Pawellek & Krivov, 2015), this $R_{\rm{BB}}$ is consistent with the
modelled inner edge of the gas. Moreover, if the CO gas in this scenario is
produced in the known planetesimal belt and observed at ${\sim}130$ au, this
backs up the previous claim that the gas must be outflowing. Thus, we may be
witnessing a short-lived phase after protoplanetary disc dispersal in which CO
ice is released and dispersed. For example, it may be that following
primordial gas dispersal a previously stable reservoir of CO ice became
susceptible to sublimation (similar to the mechanism suggested on 486958
Arrokoth by Steckloff et al., 2020), or CO was released as planetesimals were
ground down in collisions (Marino et al., 2020). If this is the case, then we
may find more examples of class III stars with rapidly dispersing gas winds.
Gas winds have not been seen towards older gaseous debris discs which may
indicate that these winds are linked to the evolutionary stage (or spectral
type) of the star, or perhaps suppressed by the build-up of other gaseous
species over many Myr.

##  6 Conclusions

By analysing the CO gas detected towards the class III star, NO Lup, we have
demonstrated that the CO has a high velocity width and is spatially resolved.
Although we showed that this can be fitted with a Keplerian Disc Model, this
requires the stellar mass of NO Lup to be implausibly high, i.e., 10 times
higher than expected. Instead, we have shown that the gas may be outflowing
with a high radial velocity, explaining the ${\sim}19$ km s-1 width and
${\sim}$130 au spatial scale. We conclude that this gas is outflowing and is
being continually replenished, and suggest possible interpretations. We note
further work to explore the nature of this detection, such as new high
resolution imaging of the continuum, measurements for the neutral CI gas, and
detailed modelling.

## Acknowledgements

We thank the anonymous reviewer for their comments which improved the quality
of this work. JBL is supported by an STFC postgraduate studentship. GMK is
supported by the Royal Society as a Royal Society University Research Fellow.
SM is supported by a Research Fellowship from Jesus College, Cambridge. GR
acknowledges support from the Netherlands Organisation for Scientific Research
(NWO, program number 016.Veni.192.233). M.T. has been supported by the UK
Science and Technology research Council (STFC), and by the European Union’s
Horizon 2020 research and innovation programme under the Marie Sklodowska-
Curie grant agreement No. 823823 (RISE DUSTBUSTERS project). JPW acknowledges
support from NSF grant AST-1907486. MK gratefully acknowledges funding by the
University of Tartu ASTRA project 2014-2020.4.01.16-0029 KOMEET. This work was
partly supported by the Deutsche Forschungs-Gemeinschaft (DFG, German Research
Foundation) - Ref no. FOR 2634/1 TE 1024/1-1. This work was partly supported
by the Italian Ministero dell’Istruzione, Università e Ricerca (MIUR) through
the grant Progetti Premiali 2012 iALMA (CUP C52I13000140001) and by the DFG
cluster of excellence Origin and Structure of the Universe (www.universe-
cluster.de).

## Data Availability

This work makes use of the following ALMA data: $ADS/$ $JAO.ALMA$
$\\#2018.1.00437.S$. ALMA is a partnership of ESO (representing its member
states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST and
ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic
of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ.
Based on observations collected at the European Southern Observatory under ESO
programme 093.C-0506(A).

## References

  * Adams et al. (1987) Adams F. C., Lada C. J., Shu F. H., 1987, ApJ, 312, 788
  * Alcalá et al. (2017) Alcalá J. M., et al., 2017, A&A, 600, A20
  * Bai (2016) Bai X.-N., 2016, ApJ, 821, 80
  * Baraffe et al. (2015) Baraffe I., Homeier D., Allard F., Chabrier G., 2015, A&A, 577, A42
  * Boccaletti et al. (2018) Boccaletti A., et al., 2018, A&A, 614, A52
  * Cieza et al. (2013) Cieza L. A., et al., 2013, ApJ, 762, 100
  * Dent et al. (2014) Dent W. R. F., et al., 2014, Science, 343, 1490
  * Ercolano & Pascucci (2017) Ercolano B., Pascucci I., 2017, Royal Society Open Science, 4, 170114
  * Foreman-Mackey et al. (2013) Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013, Publications of the Astronomical Society of the Pacific, 125, 306
  * Frasca et al. (2017) Frasca A., Biazzo K., Alcalá J. M., Manara C. F., Stelzer B., Covino E., Antoniucci S., 2017, A&A, 602, A33
  * Gaia Collaboration et al. (2018) Gaia Collaboration et al., 2018, Astronomy and Astrophysics, 616, A1
  * Hardy et al. (2015) Hardy A., et al., 2015, A&A, 583, A66
  * Haworth & Owen (2020) Haworth T. J., Owen J. E., 2020, MNRAS, 492, 5030
  * Jackson et al. (2014) Jackson A. P., Wyatt M. C., Bonsor A., Veras D., 2014, MNRAS, 440, 3757
  * Kennedy et al. (2019) Kennedy G. M., et al., 2019, Nature Astronomy, 3, 230
  * Kral et al. (2017) Kral Q., Matrà L., Wyatt M. C., Kennedy G. M., 2017, Monthly Notices of the Astronomical Society, 469, 521
  * Kral et al. (2019) Kral Q., Marino S., Wyatt M. C., Kama M., Matrà L., 2019, MNRAS, 489, 3670
  * Krautter et al. (1997) Krautter J., Wichmann R., Schmitt J. H. M. M., Alcala J. M., Neuhauser R., Terranegra L., 1997, Astronomy and Astrophysics Supplement Series, 123, 329
  * Lesur (2020) Lesur G., 2020, arXiv e-prints, p. arXiv:2007.15967
  * Lovell et al. (2020) Lovell J. B., et al., 2020, Monthly Notices of the Royal Astronomical Society
  * Manara et al. (2013) Manara C. F., et al., 2013, A&A, 551, A107
  * Marino et al. (2016) Marino S., et al., 2016, Monthly Notices of the Royal Astronomical Society, 460, 2933
  * Marino et al. (2020) Marino S., Flock M., Henning T., Kral Q., Matrà L., Wyatt M. C., 2020, MNRAS, 492, 4409
  * Matrà et al. (2015) Matrà L., Panić O., Wyatt M. C., Dent W. R. F., 2015, Monthly Notices of the Royal Astronomical Society, 447, 3936
  * Matrà et al. (2017) Matrà L., et al., 2017, Monthly Notices of the Royal Astronomical Society, 464, 1415
  * Matrà et al. (2019) Matrà L., Öberg K. I., Wilner D. J., Olofsson J., Bayo A., 2019, AJ, 157, 117
  * Modigliani et al. (2010) Modigliani A., et al., 2010, in Observatory Operations: Strategies, Processes, and Systems III. p. 773728, doi:10.1117/12.857211
  * Moór et al. (2017) Moór A., et al., 2017, ApJ, 849, 123
  * Owen & Kollmeier (2019) Owen J. E., Kollmeier J. A., 2019, MNRAS, 487, 3702
  * Pawellek & Krivov (2015) Pawellek N., Krivov A. V., 2015, Monthly Notices of the Royal Astronomical Society, 454, 3207
  * Pontoppidan et al. (2011) Pontoppidan K. M., Blake G. A., Smette A., 2011, ApJ, 733, 84
  * Siess et al. (2000) Siess L., Dufour E., Forestini M., 2000, A&A, 358, 593
  * Steckloff et al. (2020) Steckloff J. K., Lisse C. M., Safrit T. K., Bosh A. S., Lyra W., Sarid G., 2020, arXiv e-prints, p. arXiv:2007.12657
  * Visser et al. (2009) Visser R., van Dishoeck E. F., Black J. H., 2009, A&A, 503, 323
  * Williams & Cieza (2011) Williams J. P., Cieza L. A., 2011, ARA&A, 49, 67
  * Wyatt (2008) Wyatt M. C., 2008, Annual Review of Astronomy and Astrophysics, 46, 339
  * Wyatt et al. (2015) Wyatt M. C., Panić O., Kennedy G. M., Matrà L., 2015, Ap&SS, 357, 103

##  Appendix A PV Diagrams

In Table 1 we show that the two models that fit the observations have disc
position angles that are ${\sim}90^{\circ}$ different. Since our modelling
fits all channel maps in the data set (i.e., it does not fit a PV diagram) we
have a choice in how we present this. In Fig. 3 we used a slit angle of
$36^{\circ}$, i.e., parallel to direction of motion of the peak in the channel
maps, the Disc Model major axis, and the Outflow Model minor axis, which
produced a stripe (see panels 1, 2 and 4). If instead we place the PV diagram
slit PA along the major axis of the Outflow Model and the minor axis of the
Disc Model, these yield model PV diagrams that are then elliptical as can be
seen in Fig. 4, and similar to those presented in Haworth & Owen (2020). Thus,
although the choice of slit PA does not bias our modelling, a different choice
will result in a different visualisation, but not one that can distinguish
between the Keplerian Disc Model and the Outflow Model.

Figure 4: PV diagrams for the Disc Model and Outflow Model with their slits
taken at $90^{\circ}$ to those in Fig. 3 to demonstrate their elliptical
nature in this orientation. The Outflow Model is rotated slightly anti-
clockwise by the Keplerian motion, though this would be much greater if the
radial and Keplerian velocities were more similar.


# MARVIN: An Open Machine Learning Corpus and Environment for Automated
Machine Learning Primitive Annotation and Execution

Chris A. Mattmann
Jet Propulsion Laboratory
California Institute of Technology
Pasadena, CA 91109 USA
<EMAIL_ADDRESS>Sujen Shan
Jet Propulsion Laboratory
California Institute of Technology
Pasadena, CA 91109 USA
<EMAIL_ADDRESS>Brian Wilson
Jet Propulsion Laboratory
California Institute of Technology
Pasadena, CA 91109 USA
<EMAIL_ADDRESS>
###### Abstract

In this demo paper, we introduce the DARPA D3M program for automatic machine
learning (ML) and JPL’s MARVIN tool that provides an environment to locate,
annotate, and execute machine learning primitives for use in ML pipelines.
MARVIN is a web-based application and associated back-end interface written in
Python that enables composition of ML pipelines from hundreds of primitives
from the world of Scikit-Learn, Keras, DL4J, and other widely used libraries.
MARVIN allows for the creation of Docker containers that run on Kubernetes
clusters within DARPA to provide an execution environment for automated
machine learning. MARVIN currently contains over 400 datasets and challenge
problems from a wide array of ML domains including routine classification and
regression to advanced video/image classification and remote sensing.

##  1 Introduction

The world of machine learning (ML) is exploding and new ways of applying ML
for mundane, repetitive, and complex tasks is increasing: from using ML to
reply to e-mail messages automatically [1] to self-driving cars and using ML
for sustainable mobility [2] . Today, machine learning requires human capital
to perform a variety of tasks including 1) composition and aggregation of
machine learning “primitives” [3], e.g., kernels of code that clean, decide,
classify/label, cluster/organize, and rank the input dataset to produce an
output; 2) selection of ML models that have already been generated such as
from OpenML [4], ModelZoo [5], UC Irvine ML repository [6], Kaggle [7], etc.
or the generation of new models, datasets, and problems; and 3) development of
metrics to evaluate how well machine learning is working on a particular
problem.

Figure 1: D3M Government Team responsibilities.

In 2016, the United States Defense Advanced Research Projects Agency (DARPA)
incepted the Data Driven Discovery of Models (D3M) program to automate the
three core tasks of machine learning identified in the preceding paragraph.
Twenty four performers on D3M are working together in three technical areas
corresponding to each of the core ML tasks (TA1: ML primitives development;
TA2: Automated Machine Learning using Datasets, Problems, and Challenges; and
TA3: User/SME Metrics for evaluating automated machine learning). The Jet
Propulsion Laboratory, California Institute of Technology (NASA-JPL) is part
of the D3M government team, which consists of NASA-JPL, MIT Lincoln
Laboratories, and the US National Institutes of Standards and Technology
(NIST). The government team responsibilities, shown in Fig. 1, involve the
creation of a baseline open source model corpus (TA1); the curation of 10s of
thousands of ML datasets, problems, and baseline solutions (TA2) for proxy
government and even direct governmental challenges across defense, civil
agencies, and space; and a set of automated evaluation metrics to assess the
performance of the automatic machine learning systems.

Our team at NASA-JPL has over the past two years designed and implemented a
system that we call MARVIN (a reference to the character from the _Looney
Tunes_ cartoon series [8]) whose main goals are to provide an environment for
creating a performer baseline for machine learning primitive annotation,
schema, and evaluation. MARVIN is a web-based tool that is built on
ElasticSearch [9] and FacetView [10] and Kibana [11]. ElasticSearch is used to
capture JSON instances of metadata about machine learning primitives
representing the world of ML - from SKLearn [12], to DL4J [13], to Keras [14]
and so on. Primitive “annotations” capture their inputs, outputs, and
characteristics for composing the primitives into an ML pipeline for use by
the TA2 performers. MARVIN allows for search and exploration of these
primitive instances, and we have written an application programming interface
to leverage MARVIN for invocation of the associated primitives (e.g., pick
your LogisticRegression that supports particular hyperparameters, programming
languages, etc.). In addition, Marvin also provides the ability to search ML
datasets, and challenge problems for over 400 challenges at present (and soon
10s of thousands) ranging from introductory machine learning activities and
data types e.g., predict Hall of Fame decisions based on baseball player
statistics all the way to complex multi-source classification and decision
making e.g., identify objects in videos and images with high accuracy. All
primitives, datasets, and challenge problems are searchable and indexed in
ElasticSearch and browseable using FacetView and Kibana.

MARVIN also provides the ability to generate a Dockerfile and allow for
execution of ML primitives and pipelines based on the primitive type and
associated preferred data and problem domain. We have classified domains into
areas of video, text/NLP, classification, labeling, regression, etc. MARVIN
also provides base Docker images containing the curated library of ML
primitives and a base execution environment that D3M performers can extend and
use to run generated pipelines.

##  2 Data model

The Primitive “annotations” follow a [15] that allow TA2 systems to infer
metadata about the primitive, like its “algorithm type”, “primitive family”,
“hyperparameters”, etc. (JPL developed an initial schema but now it is jointly
developed by the entire D3M community – Performers and the government teams.)

Algorithm Types: The algorithm type field intends to capture the underlying
implementation of the primitive. The field currently is an enumeration of
values from a controlled vocabulary including, but not limited to, algorithms
such as ADABOOST [16], BAYESIAN LINEAR REGRESSION [17], DECISION TREE [18],
etc.

Figure 2: A screenshot of the MARVIN interface showing primitives (on the
right) and their metadata, along with search criteria on the left, e.g.,
NO_JAGGED_VALUES, etc.

Primitive Family: This field intends to capture the high level purpose of the
primitive. It is an enumeration of controlled values e.g., “CLASSIFICATION”,
“FEATURE_SELECTION”, “DATA_TRANSFORMATION”, “TIMESERIES_FORECASTING”, etc.

Hyperparameter: The hyperparameter field captures all the parameters that a
TA2 system could pass to a TA1 primitive to control its behavior. These could
include tunable hyperparameters which affect the model performance, resource
use hyperparameters that control the resource usage of that primitive, and
metafeature parameters that control which meta-features are computed by the
primitive.

Apart from the above, there are more key fields that help a TA2 system to
determine whether to use the given primitive in a pipeline for solving a
problem or not. For example, the preconditions field indicates what
preconditions must evaluate to true before using this primitive. For instance
the “NO_MISSING_VALUES” precondition would inform a TA2 system that this
primitive cannot handle missing values in the data. Similarly there could be
other preconditions like “NO_CONTINUOUS_VALUES”, “NO_JAGGED_VALUES”, etc., as
show in Fig. 2.

The effects field in the schema informs a TA2 about the set of post conditions
obtained by data processed by this primitive. Some of the effects include -
“NO_MISSING_VALUES”, this indicates that the primitive removes missing values,
“NO_CATEGORICAL_VALUES”, indicates that the primitive removes categorical
columns, etc.

##  3 Execution environment

Figure 3: Our MARVIN execution environment architecture running on top of
Kubernetes (K8s).

MARVIN leverages Docker containers and Kubernetes clusters (K8s) [19] to
provide an environment for D3M performer testing (called the “Play Space” in
Fig. 1 upper left corner). Depending on the machine learning problem type
e.g., Image/Video, NLP/Text, Classification, performers can use MARVIN to
create a Dockerfile containing a set of Primitives for a particular problem
set and then create a running Docker image on our D3M Kubernetes cluster. As
shown in Fig. 3 starting from the left side of the figure, TA1 performers can
publish a Docker image of their TA1 primitive if it requires special services
such as e.g., graph databases, or particular Python libraries which can then
be consumable and included in other performer containers. In addition, we
offer base Master containers for our TA2 systems (top-middle of Fig. 3) to
leverage e.g., NLP, Image/Video base, or a combination of all required
libraries in the program. Then, using MARVIN’s API, TA2 performers can
dynamically select, and execute their primitives. Data (shown in the middle-
left of Fig. 3) is made available via shared K8s volumes to all the containers
on the cluster.

##  4 Conclusion

We presented MARVIN, an automated system for discovering machine learning
primitives, and executing them in an automated fashion. MARVIN is a core tool
in the DARPA D3M program. At present MARVIN provides a web-based interface to
browse over hundreds of ML primitives and 400 ML datasets and challenge
problems, and over the next two years the tool will be expanded to support
thousands of ML problems and primitives for the performers in the program.

## Acknowledgement

This effort was supported in part by JPL, managed by the California Institute
of Technology on behalf of NASA, and additionally in part by the DARPA
Memex/XDATA/D3M/ASED programs and NSF award numbers ICER-1639753, PLR-1348450
and PLR-144562 funded a portion of the work. The authors would like to thank
Maziyar Boustani, Kyle Hundman, Asitang Mishra, Tom Barber and Giuseppe Totaro
at JPL who materially contributed to the design and construction of MARVIN. We
also thank our collaborators at MIT Lincoln Labs and NIST. The library of
primitives is now being further populated by the entire D3M community. We also
credit Mr. Wade Shen for his support of our work.

## References

  * [1] A. Kannan, K. Kurach, S. Ravi, T. Kaufmann, A. Tomkins, B. Miklos, G. Corrado, L. Lukács, M. Ganea, P. Young _et al._ , “Smart reply: Automated response suggestion for email,” in _Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_. ACM, 2016, pp. 955–964.
  * [2] L. D. Burns, “Sustainable mobility: a vision of our transport future,” _Nature_ , vol. 497, no. 7448, p. 181, 2013.
  * [3] T. M. Mitchell, _The discipline of machine learning_. Carnegie Mellon University, School of Computer Science, Machine Learning Department, 2006, vol. 9.
  * [4] J. Vanschoren, J. N. Van Rijn, B. Bischl, and L. Torgo, “Openml: networked science in machine learning,” _ACM SIGKDD Explorations Newsletter_ , vol. 15, no. 2, pp. 49–60, 2014.
  * [5] Y. Jia and E. Shelhamer, “Caffe model zoo,” 2015.
  * [6] A. Frank and A. Asuncion, “Uci machine learning repository [http://archive. ics. uci. edu/ml]. irvine, ca: University of california,” _School of information and computer science_ , vol. 213, 2010.
  * [7] Wikipedia, “Kaggle — Wikipedia, the free encyclopedia,” http://en.wikipedia.org/w/index.php?title=Kaggle&oldid=838984387, 2018, [Online; accessed 24-May-2018].
  * [8] J. Beck, _Looney Tunes: The Ultimate Visual Guide_. DK Publishing (Dorling Kindersley), 2003.
  * [9] C. Gormley and Z. Tong, _Elasticsearch: The Definitive Guide: A Distributed Real-Time Search and Analytics Engine_. ” O’Reilly Media, Inc.”, 2015.
  * [10] Open Knowledge Foundation, _FacetView_ , Open Knowledge Foundation, 2018. [Online]. Available: https://github.com/okfn/facetview
  * [11] Y. Gupta, _Kibana Essentials_. Packt Publishing Ltd, 2015.
  * [12] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg _et al._ , “Scikit-learn: Machine learning in python,” _Journal of machine learning research_ , vol. 12, no. Oct, pp. 2825–2830, 2011.
  * [13] S. Bahrampour, N. Ramakrishnan, L. Schott, and M. Shah, “Comparative study of caffe, neon, theano, and torch for deep learning,” 2016.
  * [14] S. Raschka, _Python machine learning_. Packt Publishing Ltd, 2015.
  * [15] D. P. A. Schema, 2017. [Online]. Available: https://gitlab.com/datadrivendiscovery/metadata/blob/master/d3m_metadata/schemas/v0/definitions.json
  * [16] Y. Freund and R. E. Schapire, “A decision-theoretic generalization of on-line learning and an application to boosting,” 1997.
  * [17] T. P. Minka, “Bayesian linear regression,” 1998.
  * [18] J. QUINLAN, “Induction of decision trees,” 1986.
  * [19] D. Bernstein, “Containers and cloud: From lxc to docker to kubernetes,” _IEEE Cloud Computing_ , vol. 1, no. 3, pp. 81–84, 2014.


Corresponding author<EMAIL_ADDRESS>
# Study on the Kinetics of Rayleigh Particle Jets Converging by Laser Beams

Z. L. Wang  [  Shanghai Institute of Applied Mathematics and Mechanics,
Shanghai Key Laboratory of Mechanics in Energy Engineering, Shanghai
University, Yanchang Road 149, Shanghai, 200072, P.R. China  Nu Zhenyu
Shanghai Institute of Applied Mathematics and Mechanics, Shanghai Key
Laboratory of Mechanics in Energy Engineering, Shanghai University, Yanchang
Road 149, Shanghai, 200072, P.R. China  Huang Kai  Shanghai Institute of
Applied Mathematics and Mechanics, Shanghai Key Laboratory of Mechanics in
Energy Engineering, Shanghai University, Yanchang Road 149, Shanghai, 200072,
P.R. China

###### Abstract

This paper discusses laser-induced flow stabilizing of Rayleigh particle jets.
Laser technology, has important applications in micro/nano-scale static
monomer particle operations, such as optical tweezers, or is used for the
passive measurement of macroscopic physical features of particle groups.
However, it is relatively rare for the laser beam to directly interfere with
the behavior of particle populations dynamically, so as to achieve the purpose
of instant group manipulations. Based on the theoretical analysis of particle
dynamics and hydrodynamic stability theory, the effects of light induced
convergence on rarified jets (as a point source emitting particle off a
nozzle) and denser jets consists of Rayleigh sized particles have been
considered. For rarified particle jet’s analysis, compared with the classical
vacuum evaporation deposition theory, we found that the laser positively
guided the movement of particles, leading their pathes into more concentrated
targets. Such convergence effect also happens in the case of denser Rayleigh
particle jets. Particle dynamics simulations and hydrodynamic stability
analysis mutually authenticated that the optical field forces suppress the
instability both of long-wave and short-wave on particle jet interfaces, and
have broad-spectrum stabilization characteristics. Therefor, diffusive
particles in vacuum evaporation can also have very good targeted aggregations
by laser.

##  I Introduction

Jets composed of particulate matter are widely present in a variety of natural
phenomena and engineering processes, such as inkjet printing, shot peening,
pulverized coal gasification, and the like. In order to study particle jets,
many methods for producing particle jets have been developed. Thoroddsen et
al. used large-diameter ($1.34cm$) high-velocity balls to impinge on a
substrate packed with $80\mu m$ diameter glass beads, demonstrating that the
particulate matters can still generate high-density jets at the micro-scale
through deformation-rebounding of the substrate Thoroddsen and Shen (2001).
Particle jets can also be formed by carrying particles through a high-speed
air stream, and a liquid film-like flow structure can also be obtained after
striking the wall surface. Shi et al. used such method to find that the
thickness of the liquid film was independent of the flow velocity and particle
size, and the evolution of the liquid film was related to the ratio of flow
diameter and particle size Shi et al. (2017).

The demands for MEMS researches have been increasing in recent years.
Researchers at Harvard University adopted $100nm$ diameter powder particles in
2013 to successfully manufacture lithium-ion micro-cells with local accuracy
of up to $1\mu m$ by using 3D printing technology Sun et al. (2013). The MIT
Microsystems Technology Laboratory developed a low-cost processing method
based on electrospray printing technology in 2015, using a solution containing
nanoparticles to produce a microsensor with only $0.03mm^{2}$ large, whose
metal wire of its internal circuit is only $10\mu m$ wide with a minimum
spacing of $50\mu m$ Taylor and Velasquez-Garcia (2015). Additive
Manufacturing technology at a micro/nano scale has also been rapidly
developed. Fuller et al. used piezoelectric ink-jet printing technology to
transport nano-particles to successfully construct a variety of micro-
electronic components, demonstrating the great potential of inkjet printing
technology in the MEMS field Fuller et al. (2002). As an alternative to inkjet
printing technology, Huang et al. used aerosol printing technology to
successfully complete the low-cost and fast printing of conductive silver
films on $A4$ paper substrates at room temperature Huang et al. (2014).

For the transportation of micro/nanoparticles in jet researches, the studies
usually focus on particles of size no smaller than $1\mu m$, and the influence
of the fluid environment is often negligible. While for nanoparticles, there
are many applications other than inkjet printing and aerosol printing. Xiao et
al. used electrospinning technology to disperse nano-metal particles into
fiber materials, which reduced the occurrence of particle agglomeration and
improved the material properties Xiao et al. (2009). Hua et al. studied the
effects of process parameters on the properties of aluminum films fabricated
on single crystal silicon by DC magnetron sputtering, and found that the films
prepared have preferred orientation for specific crystal structure particles
Hua et al. (2015). Yu et al. successfully prepared high-purity, high-
catalytic-efficiency SiC nanoparticles using heat plasma technology Yu et al.
(2017).

Because particles transportation in jets often requires various liquids,
gases, or other substances as carriers, the roles of carriers can not be
ignored, which limits their practical applications. In this paper, we attempts
to construct a theoritical basis for a direct writing method with little or
non carrier effects by using the optical field confinement effect of the
optical tweezer technique.

The study of optical tweezers originated in the $1970s$ when Ashkin in the
laboratory accidentally observed the phenomenon of laser disturbance of
colloidal particles in a solution Ashkin (1970a). Later, Ashkin and his
collaborators conducted theoretical experiments using lasers for atom/molecule
and particle capturing Ashkin (1970b, 1978). Until 1986, Ashkin successfully
used a single-beam gradient force trap for the first time in the laboratory to
achieve stable capture of colloidal particles ranging from $25nm$ to $10\mu
m$, which marked the official birth of optical tweezers Ashkin et al. (1986).
Optical tweezers technology has been widely used in microscale researches such
as chromosome separation and measurement of particle-particle interactions Ren
(2008); Wang et al. (2004).

The existing optical technology literatures focus on the researches of
manipulations of single particles. The research on the influence of a large
number of particles on population behavior is rare in literatures and patents.
While studying the vacuum evaporation coating technology of organic light-
emitting diodes (OLEDs), inspired by the phenomenon of light trapping
particles near the optical axis, this article attempts to explore the
feasibility of using optical tweezers to control the particulate flow
converging or focusing Huang et al. (2013). The force calculation method for
spherical particles in an electromagnetic field had been given in the work of
Lorenz, Mie and Debye, Lorenz-Mie theory, but their method are cumbersome to
calculate. After the birth of optical tweezers technology, in order to meet
the needs of optical tweezer designs, people had developed a variety of
methods to simplify the calculation processes. The methods applied to the
light force of the particles computing is generally differentiated according
to $\lambda/a$, the ratio of the laser wavelength to the radius of the
particles. When $\lambda/a>20$, the particles are called Rayleigh particles,
and very accurate results can be obtained using the electric dipole model
Harada and Asakura (1996). And when $\lambda/a<1/5$, particles, called Mie
particles, can be approximated by geometrical optics Ashkin (1970a). The more
generalized GLMT (Generalized Lorenz-Mie theory) method can be applied to the
calculation of spherical particle forces at various scales Gouesbet et al.
(1988).

In our previous work, the calculations of Mie particles were testified by
means of theoretical analysis and numerical simulation Huang (2013); Wang
(2016). The discussions in this article will focus on the behavior of nano-
sized Rayleigh particles under laser light beams.

(a) (a)

(b) (b)

Figure 1:  LABEL:sub@fig1a Conceptual prototype of laser induced deposition.
LABEL:sub@fig1b. Simplified computational model as particles emitting from a
point source into a Gauss laser beam.

##  II Physical model and simulation methods

###  II.1 Physical model

In this paper, the software LAMMPS(Large-scale Atomic/Molecular Massively
Parallel Simulator) was used for numerical simulation. Lammps is a molecular
dynamics software developed by the Sandia National Laboratories in the United
States. Modeling simulations of motions in different force fields and boundary
conditions such as atoms, biomolecules, particle streams, or other coarse-
grained systems can be performed. At the same time, the MPI (Message Passing
Interface) is integrated into its software architecture, simplifying the
implementation of parallel computing Plimpton (1995). And the software can be
sourced from its official website (http://lammps.sandia.gov) and authorized to
be modified when necessary.

The initial physical model can be seen in Figure 1(a) as a combination of
Knudson box (particle source) and deposition liner used in vacuum evaporation.
The laser is added to emit from the particle outflow hole, pointing to the
liner or collection device, and inducing the movement of molecular clusters or
particles in a vacuum environment.

After initially attempting to use LAMMPS to simulate the position of the
deposition liner from the inside of the evaporation chamber to the outside, it
was found that the flow simulation inside of Knudson box requires a large
amount of computational resources. And although the Knudsen number (average
molecular free path to orifice diameter ratio) is at a minimum of $0.003$,
which is close to a viscous flow state. However, the particle flow outflow and
velocity distribution are still affected by parameters such as pore length,
mean free path and pore size, and they behave more like molecular flows. And
there are too many variables that affect the simulation results Masayoshi
(2009). As the ultimate concern is the convergence effect of the laser on the
particles, taking considerations comprehensively, the model can be simplified
as combinations of a particle source and the external laser field. The
simplified model is shown in Figure 1(b). Where $r$ is the distance of a point
on the counting plane to the optical axis; $\omega_{0}$ is the waist radius of
the laser; $h$ is the coordinate of the counting plane in the $z$ direction;
$z_{0}$ is the $z$ coordinate of the waist of the laser beam.

This article will select the point source and the surface source,
respectively, to discuss the effects of laser beams on the convergence of the
Rayleigh particle jet flow.

###  II.2 Rayleigh particles affected by laser force

According to Harada’s work, the scattering force $F_{scat}$ and the gradient
force $F_{grad}$ produced by the laser beam on Rayleigh particles can be
calculated by Equation (1) Harada and Asakura (1996).

$\left\\{\begin{aligned}
{{\bf{F}}_{{\rm{scat}}}}&=\hat{z}\frac{8}{3}\pi{k^{4}}{a^{6}}\frac{{{m^{2}}-1}}{{{m^{2}}+2}}\frac{n}{c}I\left({x,y,z}\right)\\\
{{\bf{F}}_{{\rm{grad}}}}&=\frac{{2\pi
n{a^{3}}}}{c}\left({\frac{{{m^{2}}-1}}{{{m^{2}}+2}}}\right)\nabla
I(x,y,z)\end{aligned}\right.$ (1)

Where $\hat{z}$ is a unit vector in the $z$ direction, which is also the laser
propagation direction; the wave number $k=2\pi/\lambda$, where $\lambda$ is
the laser wavelength; $a$ is the radius of the particle; $m$ is the relative
refractive index of the particle to the environment; $n$ is the refractive
index of the medium in which the particles are located; $c$ is the speed of
light; and $I(x,y,z)$ is the average Poynting vector.

In this paper, the common linearly polarized Gaussian beam is selected as the
laser source. The average Poynting vector is written as Equation (2) in the
coordinate system in Fig. 1(b).

$I\left({x,y,z}\right)=\frac{{\left({\frac{{2P}}{{\pi\omega_{0}^{2}}}}\right)}}{{1+{{\left({2\frac{{z-{z_{0}}}}{{k\omega_{0}^{2}}}}\right)}^{2}}}}Exp[\frac{{-2\frac{{{x^{2}}+{y^{2}}}}{{\omega_{0}^{2}}}}}{{1+{{\left({2\frac{{z-{z_{0}}}}{{k\omega_{0}^{2}}}}\right)}^{2}}}}]$
(2)

###  II.3 Numerical aspects

The size of the simulation box set here is $100\mu m\times 100\mu m\times
81\mu m$. Different particle sources are placed at the origin of the
coordinates of Figure 1(b), emitting particles in the positive $z$-axis
direction.

Rayleigh particles here were considered equivalent to ideal spheres in our
simulations, without introducing extra complexity of their specific
morphology. Particles move in an almost vacuum environment. Therefore, the
main considerations are the interaction of particles with each other, and also
the interactions between particle and laser beam.

According to the setting scale of the model, this paper deals with a mesoscale
problem. The commonly used particle methods on such scales include multi-
particle collision dynamics (MPC), dissipative particle dynamics (DPD), and
direct simulation Monte Carlo (DSMC), etc. Noguchi et al. (2007). The DPD
method is often used to study the movement of macromolecules in fluids. By
applying coarse-grained modeling of fluids and molecules to the applied
potential obtained by statistical methods, the amount of calculations can be
reduced to speed up the simulation process. The coarse-grained process makes
the potential of the DPD method naturally include a part that reflects the
viscous and thermal motion of the fluid environment Groot and Warren (1997).
MPC and DSMC methods are similar to DPD, have force terms counting the
influence of the fluid environment. However, the objects that these methods
target are different from the goals that this article is trying to simulate.

Referring to the work of Israelachvili, Israelachvili (1985), the attractive
force between nano-sized homogeneous spherical particles decays rapidly within
a surface spacing much smaller than the diameter, and it attracts each other
during extremely close distances and collisions. It is difficult to calculate
due to the deformation of the particles, and the overall repulsive force
manifests itself in the elasto-plastic deformation of the particles when they
collide. At present, no description has been found for the particle
attractiveness that fits this problem. According to Hamaker’s work, the mutual
attraction between the particles in the simulation is simulated, and the
attractive force can approach the laser force when the particle surface
spacing is smaller than the diameter Hamaker (1937). For the repulsion between
nanoparticles, considering the particle flow study, due to the observation of
particles larger than $1\mu m$, the mutual attraction between particles due to
Van der Waals forces can be ignored. The particle collision process is
calculated according to the extrusion deformation, in which the plastic
deformation occurs. The impact of mass changes, etc., is reflected in
collisional dissipation Goldhirsch (2003). Therefore, in this paper, the
impact of particle flow collision developed by Silbert and Brilliantov was
used in the simulation to simulate the repulsion between nanoparticles Silbert
et al. (2002); Brilliantov et al. (1996), but the weak mutual attraction
between nanoparticles was not simulated. In this paper, the simulation of the
interaction between nanoparticles is not perfect, and further work is needed
for discussion.

At this time, the interaction between the spherical particles i,j is
calculated according to the collisional effect of the particle stream, ie,
equation (3).

${{\bf{F}}_{{\rm{ij}}}}=\left({{k_{n}}\delta{n_{ij}}-{m_{eff}}{\gamma_{n}}{v_{n}}}\right)-({k_{t}}{\rm{\Delta}}{s_{t}}+{m_{eff}}{\gamma_{s}}{v_{s}})$
(3)

Where $k_{n}$ is the normal elastic coefficient between the particles (ball
center line direction); $\delta n_{ij}$ is the variation of the center
distance of the particle; $m_{eff}=m_{i}m_{j}/(m_{i}+m_{j})$ is the effective
mass at the time of particle collision; $\gamma_{n}$ is the normal damping
Coefficient; $v_{n}$ is the component of the relative speed of two particles
in the normal direction; $k_{t}$ is the tangential elasticity coefficient;
$\Delta s_{t}$ is the tangential displacement vector of two balls;
$\gamma_{s}$ tangential damping coefficient; $v_{s}$ is the tangential
component of relative velocity. Among them, the elastic coefficient and the
damping coefficient need to be manually set. Since the effects of collisional
dissipation and particle rotation are not considered in this paper,
$k_{t}=0,k_{n}=0.005nN/\mu m$ are set in the simulation, and the damping
coefficient calculation is turned off.

In the simulation, the particle motion trajectory is obtained by integrating
time (4) and equation (5), that is, the Velocity-Verlet algorithm Frenkel and
Smit (1997).

${\bf{r}}\left({t+{\rm{\Delta}}t}\right)={\bf{r}}\left(t\right)+{\bf{v}}(t){\rm{\Delta}}t+\frac{{{\bf{F}}\left(t\right)}}{{2m}}{\rm{\Delta}}{t^{2}}$
(4)

${\bf{v}}\left({t+{\rm{\Delta}}t}\right)={\bf{v}}\left(t\right)+\frac{{{\bf{F}}\left({t+{\rm{\Delta}}t}\right)+{\bf{F}}\left(t\right)}}{{2m}}{\rm{\Delta}}t$
(5)

Where $r$ is the particle position vector, $v$ is the particle velocity
vector, $F$ is the particle force, $t$ is the current time, and $\Delta t$ is
the time step.

In this paper, glycerol, a common test substance, was selected as the source
of particle parameters for numerical simulation, and the environmental
assumption of vacuum evaporation was used. The particle density
$\rho=1.261×10^{3}kg/m^{3}$, the relative refractive index $m=1.4746$, and the
particle radius a is $0.03\mu m$. Ambient refractive index $n=1$(vacuum).

##  III Different particle source simulations and results

###  III.1 Point source simulation

This article first uses the point source model to begin research, preliminary
verification procedures and the feasibility of laser beam collection
particles. The initial velocity of the particles needs to satisfy a certain
velocity distribution. For this reason, reference is made herein to the
description of various evaporation sources in the vacuum evaporation
technique. In many vacuum evaporation processes, particles do not almost
collide during the flight except at the outlet of the evaporation source
Ohring (2001). This article chooses a point source that is easy to get the
velocity distribution as a starting point. The particle velocity probability
density distribution from the point particle source is isotropic. Based on
this property and the basic principles of statistical mechanics, it is
possible to give a system of equations that satisfies the probability density
function $f(v)$ of the velocity vector v in the polar coordinate space when
the particle mean initial kinetic energy $\overline{\epsilon}$ is constant. 6)
Further, the probability density distribution of the initial velocity of the
particle originated from the point source is derived as formula (7), and the
probability density distribution satisfying each velocity component is given
by (8) Zhao and Luo (2005).

$\left\\{\begin{aligned}
\mathop{\smallint}\limits_{0}^{2\pi}{\rm{d}}\phi\mathop{\smallint}\limits_{0}^{\pi}\sin\theta{\rm{d}}\theta\mathop{\smallint}\limits_{0}^{\infty}{v^{2}}{\rm{d}}vf(v,\theta,\phi)&=1\\\
\frac{1}{2}m\mathop{\smallint}\limits_{0}^{2\pi}{\rm{d}}\phi\mathop{\smallint}\limits_{0}^{\pi}\sin\theta{\rm{d}}\theta\mathop{\smallint}\limits_{0}^{\infty}{v^{4}}{\rm{d}}vf(v,\theta,\phi)&=\overline{\epsilon}\end{aligned}\right.$
(6)
${\rm{f}}\left(v\right)={\left({\frac{{3m}}{{4\pi\overline{\epsilon}}}}\right)^{\frac{3}{2}}}Exp\left[{-\frac{{3m{v^{2}}}}{{4\overline{\epsilon}}}}\right]$
(7)
${{\rm{f}}_{i}}\left({{v_{i}}}\right)={\left({\frac{{3m}}{{4\pi\overline{\epsilon}}}}\right)^{\frac{1}{2}}}Exp\left[{-\frac{{3mv_{i}^{2}}}{{4\overline{\epsilon}}}}\right],(i=x,y,z).$
(8)

With reference to Oring’s work on vacuum deposition deposition thickness, it
can be obtained on the counting plane shown in Fig. 1(b). When no laser beam
is collected, the density of the number of particles passing through the plane
at the distance $r$ from the optical axis is satisfied. Dimensionless
Distributed (9) Ohring (2001).

$\frac{{\rm{N}}}{{{{\rm{N}}_{0}}}}=\frac{1}{{{{\left({1+{{\left({\frac{r}{h}}\right)}^{2}}}\right)}^{\frac{3}{2}}}}}.$
(9)

Where $N$ is the number of crossings of the particle from the optical axis
$r$; the number of crossings of the particle at the $N_{0}$ optical axis.

$300,000$ particles were randomly injected into the calculation area from the
origin of the coordinates in a time range of $0.42s$. The particle velocity is
assigned according to the velocity distribution (8) when entering the
simulation zone. Although the point source emission particles are
omnidirectional, only the particles moving upward hemisphere as shown in Fig.
1(b) are simulated here. In addition, the lower limit of the emission distance
should be set during the simulation to ensure that the newly injected
particles will not collide with the already-launched particles at the time of
injection, so that the energy cannot be conserved during the collision
calculation due to the excessive distance. In order to reduce the number of
simulation steps, the allowable time step $\Delta t$ is automatically adjusted
in the course of the simulation according to the maximum velocity of the
particles and the defined single-step maximum movement distance in the
interval of $0.0005\mu s$ to $0.05\mu s$.

Point particle source simulation laser wavelength $\lambda=632.8nm$, beam
waist position $z_{0}=5\mu m$, waist waist radius $\omega_{0}=5\mu m$. The
average initial kinetic energy of point particle source simulation is
$\overline{\epsilon}=0.0388eV$.

The cross-over density distribution at $h=10\mu m$ under the action of
different power lasers is shown in FIG. 2 . Fig. 2(a) is the original
calculation result of crossing density, and Fig. 2(b) is the result of
dimensionless reference formula (9).

(a) (a)

(b) (b)

Figure 2:  Comparison of the distribution and theoretical distribution of the
crossing density at the counting plane $h=10$ under different power
influences. LABEL:sub@fig2a. Density distribution at different powers.
LABEL:sub@fig2b. Non-dimensionalized distribution.

From FIG. 2 , it can be seen that as the laser power increases, the
distribution of the particles is concentrated near the optical axis of the
laser ($r=0$). In the absence of a laser effect ($P=0.00W$), Fig. 2(b) shows
that the simulation results are returned to the dimensionless distribution
form that should be followed in the case of no bundles calculated by equation
(9). Therefore, for the point source, the laser can effectively converge the
trajectory of the particles in the space under certain parameters, improve the
distribution in the specified plane, and at the same time, there is no
abnormality in the movement of the particles after a few collisions in the
simulation process. This program can simulate the collision and movement of
particles normally.

###  III.2 Theoretical analysis and simulation of light converging
nanoparticle jets

####  III.2.1 Light Convergence Particle Jet Stability

Lun gave the control equation (10) of the hydrodynamic form of the particle
flow in 2006 Lun (2006).

$\left\\{\begin{aligned}
\frac{{\rm{D}}}{{{\rm{D}}t}}\rho&=-\rho\nabla\bullet{\bf{V}}\\\
\rho\frac{{{\rm{D}}{\bf{V}}}}{{{\rm{D}}t}}&=\rho{\bf{F}}-\nabla\bullet{\bf{P}}\\\
\frac{3}{2}\rho\frac{{{\rm{DT}}}}{{{\rm{D}}t}}&=-\nabla\bullet Q-P:\nabla
V-\Gamma\end{aligned}\right.$ (10)

Where $\rho$ is the local average mass density, $T$ is the particle
temperature, $Q$ is the heat flow vector of the particle stream, $P$ is the
pressure tensor of the particle stream, $\Gamma$ is the energy dissipation
rate per unit volume, and $F$ is the unit mass of the acceptor force. The
particle temperature $T$ in the particle stream does not refer to the
thermodynamic temperature of the material and is the average measure of the
pulsation velocity. For specific calculations, see equation (11).

$\left\\{\begin{aligned} n\left({r,t}\right)&=\mathop{\smallint}\nolimits
f\left({v,r,t}\right)dv\\\
V\left({r,t}\right)&=\frac{1}{{n\left({r,t}\right)}}\mathop{\smallint}\nolimits
vf\left({v,r,t}\right)dv\\\
{\rm{T}}\left({r,t}\right)&=\frac{1}{{n\left({r,t}\right)}}\mathop{\smallint}\nolimits{\left({v-V}\right)^{2}}f\left({v,r,t}\right)dv\end{aligned}\right.$
(11)

Where $f(v,r,t)$ is the particle velocity distribution function, $n(r,t)$ is
the number density of particles, and $V(r,t)$ is the velocity field
distribution.

The $Q$, $P$, and $\Gamma$ of the Lun structure are all functions of the
volume fraction $\nu=\rho/\rho_{p}$ (ratio of the mass density of the jet to
the particle density) and the particle temperature. The specific form is
relatively complex and can be found in the corresponding literature Lun
(2006).

In previous work, the action of the laser was introduced by the volume force
term $F$, ie equation (12).

$\left\\{\begin{aligned} {\bf{F}}&=-\nabla\bullet\mathord{\buildrel{\lower
3.0pt\hbox{$\scriptscriptstyle\leftrightarrow$}}\over{\cal
T}}-\frac{{\partial\mathord{\buildrel{\lower
3.0pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over{g}}}}{{\partial t}}\\\
\mathord{\buildrel{\lower
3.0pt\hbox{$\scriptscriptstyle\leftrightarrow$}}\over{\cal
T}}&=-\left({{\overline{\epsilon}_{0}}EE+\frac{1}{{{{\rm{\mu}}_{0}}}}BB}\right)+\frac{1}{2}({\overline{\epsilon}_{0}}{{\rm{E}}^{2}}+\frac{1}{{{{\rm{\mu}}_{0}}}}{{\rm{B}}^{2}})\mathord{\buildrel{\lower
3.0pt\hbox{$\scriptscriptstyle\leftrightarrow$}}\over{j}}\\\
\mathord{\buildrel{\lower
3.0pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over{g}}&={\overline{\epsilon}_{0}}E\times
B\end{aligned}\right.$ (12)

Where $\mathord{\buildrel{\lower
3.0pt\hbox{$\scriptscriptstyle\leftrightarrow$}}\over{\cal T}}$ is the
momentum flow density of the electromagnetic field and
$\mathord{\buildrel{\lower
3.0pt\hbox{$\scriptscriptstyle\rightharpoonup$}}\over{g}}$ is the momentum
density of the electromagnetic field Huang (2013). $E$ and $B$ are the
electric field and magnetic field vectors, respectively.
$\mathord{\buildrel{\lower
3.0pt\hbox{$\scriptscriptstyle\leftrightarrow$}}\over{j}}$ is a unit tensor.

The following basic flow was used for stability analysis, where density
$\bar{\rho}={\rho_{0}}=const$, pressure distribution $\bar{p}=\bar{p}(r)$,
temperature distribution $\bar{T}=\bar{T}(r)$, outlet jet radius
$\bar{h}={R_{0}}$, outlet flow rate
${\bf{\bar{V}}}=(0,0,{\bar{V}_{\rm{z}}}=U=const)$, the free surface boundary
radius ${R_{0}}$ at the outlet, the outlet volume fraction
$\bar{\nu}={\nu_{0}}={\rho_{0}}/{\rho_{p}}=const$, and ${\rho_{p}}$ is the
particle density. Then superimpose the small perturbation on the elementary
flow and substitute the non-dimensionalized particle flow control equation to
obtain the perturbation control equation (13).

$\displaystyle\left\\{\begin{aligned}
\frac{{\partial\hat{\rho}}}{{\partial\hat{t}}}+\frac{{\partial\hat{\rho}}}{{\partial\hat{z}}}&=-\hat{\nabla}\bullet{\bf{\hat{V}}}\\\
\frac{{\partial{{\hat{V}}_{z}}}}{{\partial\hat{t}}}+\frac{{\partial{{\hat{V}}_{z}}}}{{\partial\hat{z}}}&={\hat{a}_{F}}{\hat{\bar{F}}_{z}}\hat{\rho}-\frac{{\partial\left({{{\hat{a}}_{p}}\hat{\rho}+{{\hat{b}}_{p}}\hat{T}}\right)}}{{\partial\hat{z}}}+\frac{2}{{{\rm{Re}}}}\frac{\partial}{{\partial\hat{z}}}\frac{{\partial{{\hat{V}}_{z}}}}{{\partial\hat{z}}}+\frac{1}{{{\rm{Re}}}}\frac{1}{{\hat{r}}}\frac{\partial}{{\partial\hat{r}}}\hat{r}(\frac{{\partial{{\hat{V}}_{z}}}}{{\partial\hat{t}}}+\frac{{\partial{{\hat{V}}_{r}}}}{{\partial\hat{z}}})+\frac{1}{{{\rm{R}}{{\rm{e}}_{\rm{\lambda}}}}}\frac{\partial}{{\partial\hat{z}}}\left({\hat{\nabla}\bullet{\bf{\hat{V}}}}\right)\\\
\frac{{\partial{{\hat{V}}_{r}}}}{{\partial\hat{t}}}+\frac{{\partial{{\hat{V}}_{r}}}}{{\partial\hat{z}}}&={\hat{a}_{F}}{\hat{\bar{F}}_{r}}\hat{\rho}-\frac{{\partial\left({{{\hat{a}}_{p}}\hat{\rho}+{{\hat{b}}_{p}}\hat{T}}\right)}}{{\partial\hat{r}}}+\frac{1}{{{\rm{Re}}}}\frac{\partial}{{\partial\hat{z}}}\left({\frac{{\partial{{\hat{V}}_{z}}}}{{\partial\hat{t}}}+\frac{{\partial{{\hat{V}}_{r}}}}{{\partial\hat{z}}}}\right)+\frac{2}{{{\rm{Re}}}}\frac{1}{{\hat{r}}}\frac{\partial}{{\partial\hat{r}}}(\hat{r}\frac{{\partial{{\hat{V}}_{r}}}}{{\partial\hat{r}}})-\frac{1}{{{\rm{Re}}}}\frac{{{{\hat{V}}_{r}}}}{{{{\hat{r}}^{2}}}}+\frac{1}{{{\rm{R}}{{\rm{e}}_{\rm{\lambda}}}}}\frac{\partial}{{\partial\hat{r}}}\left({\hat{\nabla}\bullet{\bf{\hat{V}}}}\right)\\\
\frac{{\partial\hat{T}}}{{\partial\hat{t}}}+\frac{{\partial\hat{T}}}{{\partial\hat{z}}}&=\frac{1}{{{\rm{RePr}}}}\left({\frac{\partial}{{\partial\hat{z}}}\left({\frac{{\partial\hat{T}}}{{\partial\hat{z}}}}\right)+\frac{1}{{\hat{r}}}\frac{\partial}{{\partial\hat{r}}}(\hat{r}\frac{{\partial\hat{T}}}{{\partial\hat{r}}})}\right)-{\sigma_{T}}\hat{\bar{T}}\left({\hat{\nabla}\bullet{\bf{\hat{V}}}}\right)+\left({{{\hat{a}}_{\rm{\Gamma}}}\hat{\rho}+{{\hat{b}}_{\rm{\Gamma}}}\hat{T}}\right)\end{aligned}\right.$
(13)

The Reynolds number in the perturbation control equation
${\rm{Re}}=(U{R_{0}}{\rho_{0}})/{\rm{\bar{\varsigma}}}$, the additional
Reynolds number
${\rm{R}}{{\rm{e}}_{\rm{\lambda}}}=(U{R_{0}}{\rho_{0}})/\bar{\lambda}$, the
electro-Eula number
${{\rm{E}}_{\rm{u}}}=(0.25{\epsilon_{0}}E_{0}^{2})/({\rho_{0}}{U^{2}})$, The
Prandtl number ${\rm{Pr}}=1.5{\rm{\bar{\varsigma}}}/\bar{K}$, where
$\epsilon_{0}$ is the dielectric constant of the environment and $E_{0}$ is
the maximum electric field field strength in the laser light field. Average
particle viscosities ${\rm{\bar{\varsigma}}}$ and $\bar{\lambda}$, average
thermal conductivity $\bar{K}$, and other coefficients ${\hat{a}_{p}}$,
${\hat{b}_{p}}$, ${\sigma_{T}}$, ${\hat{a}_{\rm{\Gamma}}}$, ${\hat{a}_{F}}$,
are more complex and can be found in the literature Huang (2013).

The solution (14) of the elementary stream temperature field
${\rm{\hat{\bar{T}}}}\left({\hat{r}}\right)$ and the radial force
${\hat{\bar{F}}_{r}}\left({\hat{r}}\right)$ can be obtained from equation
(13), where the form of ${\hat{a}_{t}}$, ${\hat{b}_{t}}$,
$\bar{\Theta}\left(r\right)$) and ${\Re^{\left(1\right)}}(r)$ can be found in
the corresponding literature Huang (2013).

$\left.\begin{aligned}
{\rm\hat{\bar{T}}}\left({\hat{r}}\right)&=\frac{{{{\hat{b}}_{t}}\hat{\bar{F}}{{\left(1\right)}^{2}}}}{{\hat{r}\hat{\bar{F}}{{\left({\hat{r}}\right)}^{2}}\left({-1+{{\hat{b}}_{t}}\hat{\bar{F}}{{\left(1\right)}^{2}}\mathop{\smallint}\nolimits_{1}^{\hat{r}}-\frac{{2{{\hat{a}}_{t}}}}{{\hat{\bar{F}}{{\left(t\right)}^{3}}{t^{2}}}}{\rm{d}}t}\right)}}\\\
{\hat{\bar{F}}_{r}}\left({\hat{r}}\right)&={E_{u}}\bar{\epsilon}{\hat{a}_{0}}\bar{\Theta}\left(r\right)[\frac{\partial}{{\partial\hat{r}}}[{\Re^{\left(1\right)}}{\left({\rm{r}}\right)^{2}}]]\end{aligned}\right.$
(14)

The following stream of particles was then constructed for numerical analysis.
Under a vacuum environment, the jet of particles exits from a nozzle with a
radius ${R_{0}}=3mm$ at a rate of $U=100m/s$. The outlet volume fraction of
the jet is $\nu=0.1$ and the maximum volume fraction $\nu_{m}=0.55$. The
diameter of the particles was $100nm$, the density was $1.5g/c{m^{3}}$, and
the refractive index was $1.3$. Near-elastic collisions occur between the
particles, and the elastic recovery coefficient $\eta=0.95$. The particle
temperature of the jet is $10K$. Since the focus is on the influence of the
laser on the particle jet, the inviscidity assumption
(${\rm{Re}},{\rm{R}}{{\rm{e}}_{\rm{\lambda}}}\to\infty$) is introduced to
simplify the subsequent analysis.

Referring to Drazin’s work Drazin and Reid (1984), use the regular modal
method to write the disturbance variables
${\rm{\hat{\rho}}},{\hat{V}_{z}},{\hat{V}_{r}},{\rm{\hat{T}}}$ in the form of
equation (15). Where ${\rm{\Omega}}$ is the time growth factor and $k$ is the
axial wave number. The flow is stable when ${\rm{\Omega}<0}$, and the flow is
neutrally stable when ${\rm{\Omega}=0}$, and the flow will destabilize with
time when ${\rm{\Omega}>0}$.

$\left\\{{\hat{\rho},{{\hat{V}}_{z}},{{\hat{V}}_{r}},\hat{T}}\right\\}=\\{\hat{\rho}(r),{\hat{V}_{z}}(r),{\hat{V}_{r}}(r),\hat{T}(r)\\}exp({\rm{\Omega}}\hat{t})exp(ik\hat{z})$
(15)

It is substituted into the equation set (13) and simplified using the inviscid
assumption to obtain a second-order ordinary differential equation.

$\frac{{{\partial^{2}}}}{{\partial{{\hat{r}}^{2}}}}\left({\hat{r}{{\hat{V}}_{r}}}\right)+{A_{1}}\left(r\right)\frac{\partial}{{\partial\hat{r}}}\left({\hat{r}{{\hat{V}}_{r}}}\right)+{A_{0}}\left(r\right)\left({\hat{r}{{\hat{V}}_{r}}}\right)=0$
(16)

where, ${A_{0}}\left(r\right)=({\rm{\Omega}}+ik)/{\rm{\Psi}}(r)$,
${A_{1}}\left(r\right)=-[{{\hat{a}}_{F}}{{\hat{\bar{F}}}_{r}}{B_{\rho}}\left({\rm{r}}\right)-\frac{{\partial{\rm{\Psi}}\left(r\right)}}{{\partial\hat{r}}}+\frac{1}{{\hat{r}}}{\rm{\Psi}}\left(r\right)]/{\rm{\Psi}}\left(r\right)$,
${\rm{\Psi}}\left(r\right)=\left({{{\hat{a}}_{p}}+{{\hat{b}}_{p}}{B_{T}}\left(r\right)}\right){B_{\rho}}(r)$,
${B_{\rho}}\left(r\right)=-{\rm{\Omega}}+ik+ik{B_{{V_{z}}}}{\left(r\right)^{-1}}$,
${B_{{V_{z}}}}\left(r\right)=\left({{{\hat{a}}_{F}}{{\hat{\bar{F}}}_{r}}-ik{{\hat{a}}_{p}}-ik{{\hat{b}}_{p}}{B_{T}}\left(r\right)}\right)/({\rm{\Omega}}+ik)$,
and
${B_{T}}\left(r\right)=\left({{\sigma_{T}}\hat{\bar{T}}\left({{\rm{\Omega}}+ik}\right)-{{\hat{a}}_{\rm{\Gamma}}}}\right)/({\rm{\Omega}}+ik+{{\hat{b}}_{\rm{\Gamma}}})$.

(a) (a)

(b) (b)

Figure 3: Hydrodynamic stability analysis of light forces on particle jets.
LABEL:sub@fig3 Normalized basic state radial force $\bar{F}_{r}$ and particle
temperature $\bar{T}$. LABEL:sub@fig4 Effect of electro-Euler number $Eu$ on
flow stability.

Then, the new equation (16) is again transformed into a linear system of
equations (17) using a pseudo-spectral configuration method for numerical
solution. The system of equations is a homogeneous determinant, and the
existence of a nontrivial solution is that the determinant of the matrix is
zero.

$\Xi\left({{\rm{\Omega,}}k;OtherPars}\right)=0$ (17)

The equation set (17) was solved using the Müller method Moler and Stewart
(1973). Figure 3 reflects the calculated trends of radial force and particle
temperature. In the figure, the radial force always points to the optical axis
and increases with the increase of the radius, so that the particles converge
toward the optical axis and the jets tend to converge. The analytical
expression of $\bar{F}_{r}$ also shows that the laser power is increased and
the binding force is increased. The distribution of particle temperature shows
that the particle’s activity in the vicinity of the boundary decreases rapidly
while the internal change is relatively small. Fig. 4 is the relationship
between the axial wave number $k$ and the growth factor
${{\rm{\Omega}}_{\rm{r}}}$ under different electro Euler numbers $E_{u}$. The
electron Euler number $E_{u}$ decreases in the direction of the arrow in the
figure, and the growth factor increases. From the situation change in the
figure, it can be seen that, as the Euclidean number $E_{u}$ decreases as
${E_{u}}\leq 0.01$, the growth factor ${{\rm{\Omega}}_{\rm{r}}}>0$, the
analysis system according to the regular modal method will be unstable, and
the stable jet cannot be maintained; and when the Eulerian When the number
${E_{u}}\geq 0.5$, the system remains stable. The most direct way to increase
the Eulerian number $E_{u}$ is to increase the laser power, which is
consistent with the concept of using the laser to stabilize the jet. Figure 4
also shows that the laser suppresses both the long-wave and short-wave
instability of the jet. This is different from the typical surface force in
the common stability theory—the selectivity of the surface tension on the
disturbance wave. Usually the surface tension has a destabilizing effect on
the long-wave surface of the jet, and it has an inhibitory effect on the
short-wave.

####  III.2.2 Particle Dynamics Simulation of Nanoparticle Jet

In the real scientific research or production process, the particle source may
be nozzles, helium and other components, and have a certain export area. For
the velocity distribution of particles at the outlet of the nozzle in gas jets
or jets at the macroscopic scale, although Anderson, Nathanson, et al. have
done some theoretical and experimental measurement work, but mostly focused on
the axial direction of the nozzle combined with one-dimensional theoretical
results. Discussion, lacking a more complete theoretical description Nathanson
(2004); Anderson and Fenn (2004). In the simulation process of nozzles on the
micro-nano scale, it was found that the adjustment of the geometric parameters
would cause complex changes in the outlet flow velocity distribution. Based on
the above considerations, this paper approximates the divergence of the
surface source by superimposing a disturbance velocity distribution when the
particle stream enters the simulation area. In the following simulation
results, the perturbation velocity uses the same probability density
distribution function as the aforementioned point particle source. The
perturbation velocity distribution parameter is set to
$\overline{\epsilon}=0.004eV$. The laser beam waist position $z_{0}$ was
adjusted to $20\mu m$, and the other parameters were the same as those of the
spot particle source simulation.

Figure 4: screenshots of the spatial distribution of particle jets at
different powers. (a). $P=0.0W$. (b). $P=0.5W$. (c). $P=1.0W$.

In order to reduce the amount of calculations during simulation, this paper
refers to the particle injection method by Moseler et al. Moseler and Landman
(2000). First, simulate a group of $300,000$ particles randomly moving in a
cylindrical tube with a diameter of $1\mu m\times 20\mu m$ as the liquid
column, and fix the relative positions of the particles. The liquid column was
injected at a fixed speed of $v_{z}=0.1m/s$, the relative positions of the
particles entering the simulation region were released and the disturbance
velocity was superimposed on the $v_{z}$ basis.

Figure 5: Trajectories of $50$ adjacently numbered particles at $P=1.0W$.
Figure 6: Non-dimensional cumulative distributions of particle jeta at
$h=80\mu m$ counting plane at different powers. Figure 7: Dispersive
distribution variations on $h=80\mu m$ counting plane according to changes of
electro-Euler number.

Fig. 5 is the X-directional projection of the spatial distribution of a
particle flow at a time with different power ($P=0.0W,0.5W,1.0W$) Laser band.
The dotted line in the figure is shown as the boundary of the laser beam. It
can be seen that the spatial distribution density of particles is obviously
changed. With the increase of power, the particle density in fig. 5 is
gradually reduced. This is consistent with the conclusion of the Huang Kai,
and increasing the laser power (electro-Euler number) can make the system to
form a stable jet. At the same time, there is no obvious change in the
diameter of the jet in the simulated region, that is, the instability
phenomenon in the traditional jet theory. This is consistent with the
conclusion that the Huang Kai laser will suppress the long-wave and the short
wave instability. Fig. 6 is the flight trajectory of $50$ adjacent numbered
particles in $P=1.0W$. In this paper, the collision is judged by calculating
the angle change of the flight trajectory in two adjacent time steps of the
particle after the simulation. In the simulation of Fig. 5 (c) $P=1.0W$, the
average collision number of particles in the simulated domain is $12.2$ times,
the $10^{\circ}$ time is 5.66 times, and the average collision frequency is
about $15.3\times{10^{3}}Hz$ and $7.08\times{10^{3}}Hz$ respectively
$5^{\circ}$ (The collision frequency of the molecules in the air is
${10^{9}}{10^{10}}{\rm{Hz}}$ level, i.e. $10^{2}$ order of magnitude for
velocity and $10^{-7}$ order of magnitude for molecular free path). The
particles that fly out of the beam range are about $2.36\%$ of the total
grains.

On the $h=80\mu m$ plane, statistics of cumulative distribution densities at
different powers are obtained in FIG. 7 , and it is found that the
distribution of particles in the beam range under the simulated power is not
significantly different. Use equation (18) to define the mean square radius
$\bar{r}$ to represent the degree of dispersion of the particle position in
the entire plane and observe the laser beam convergence effect. It can be seen
from the curve fitted in the figure that with the linear growth of the
electro-Euler number, the dispersion degree tends to decrease similarly to a
negative index.

##  IV Conclusion

In this paper, the dynamic simulation of the superposition of the laser light
field by the nanoparticle flow generated by the point source and the jet is
realized based on the Lammps software. The stability of the nanoparticle jet
under laser irradiation is analyzed by theoretical analysis. Comparing the
point source simulation results with the theoretical distribution of classical
vacuum evaporation deposition, the enhancement of the laser power $P$ can
increase the concentration of the particle distribution. It is initially
proved that the laser can effectively influence the particle motion
trajectory, thereby improving the deposition and distribution of the nano-
scale particles.

The theoretical analysis of the stability of photonic nanoparticle jets shows
that increasing the electro Euler number $Eu$ can improve the stability of the
particle jet to form a stable jet. At the same time, it is found that the
laser suppresses the long-wave and short-wave stability of the particle jet at
the same time. The simulation results of particle jets show that increasing
electro Euler number $Eu$ by increasing the laser power $P$ can indeed
constrain the particle formation jet more effectively, and the mean square
radius decreases with the increase of the laser power, indicating that the
concentration of particles also follows. The increase in electro-Euler number
$Eu$ has increased, demonstrating a trend of stability that is consistent with
theoretical predictions.

###### Acknowledgements.

This work was supported by National Natural Science Foundation of China
(11172163), and also partially by the Shanghai Natural Science Foundation
(15ZR1416400) and the National Key Research and Development Program
(2017YFB0404503).

## References

  * Thoroddsen and Shen (2001) S. T. Thoroddsen and A. Q. Shen, Physics of Fluids 13, 4 (2001).
  * Shi et al. (2017) Z. H. Shi, W. F. Li, W. W. Qian, H. F. Liu, and F. C. Wang, Chemical Engineering Science 162, 1 (2017).
  * Sun et al. (2013) K. Sun, T. S. Wei, B. Y. Ahn, J. Y. Seo, S. J. Dillon, and J. A. Lewis, Advanced Materials 25, 4539 (2013).
  * Taylor and Velasquez-Garcia (2015) A. P. Taylor and L. F. Velasquez-Garcia, Nanotechnology 26 (2015).
  * Fuller et al. (2002) S. B. Fuller, E. J. Wilhelm, and J. M. Jacobson, Journal of Microelectromechanical Systems 11, 54 (2002).
  * Huang et al. (2014) Q. J. Huang, W. F. Shen, Q. S. Xu, R. Q. Tan, and W. J. Song, Materials Letters 123, 124 (2014).
  * Xiao et al. (2009) S. L. Xiao, M. W. Shen, R. Guo, S. Y. Wang, and X. Y. Shi, Journal of Physical Chemistry C 113, 18062 (2009).
  * Hua et al. (2015) Y. Q. Hua, A. C. Zhu, R. F. Chen, and L. Q. Guo, Journal of Functional Materials 46, 4071 (2015).
  * Yu et al. (2017) J. Y. Yu, H. Huang, J. Gao, L. Zhou, S. Gao, X. L. Dong, and Y. Quan, Journal of Inorganic Materials 32, 351 (2017).
  * Ashkin (1970a) A. Ashkin, Phys.rev.lett 24, 156 (1970a).
  * Ashkin (1970b) A. Ashkin, Physical Review Letters 25, 1321 (1970b).
  * Ashkin (1978) A. Ashkin, Physical Review Letters 40, 729 (1978).
  * Ashkin et al. (1986) A. Ashkin, J. M. Dziedzic, J. E. Bjorkholm, and S. Chu, Optics Letters 11, 288 (1986).
  * Ren (2008) H. L. Ren, Chinese Journal of Lasers 35, 151 (2008).
  * Wang et al. (2004) H. W. Wang, X. H. Liu, Y. M. Li, B. Han, L. R. Lou, and K. J. Wang, ACTA BIOPHYSICA SINICA 20, 50 (2004).
  * Huang et al. (2013) K. Huang, Z. L. Wang, J. S. Zhou, Z. W.and Zhang, and J. H. Zhang, Journal of Shanghai University(Natural Science Edition) 19, 598 (2013).
  * Harada and Asakura (1996) Y. Harada and T. Asakura, Optics Communications 124, 529 (1996).
  * Gouesbet et al. (1988) G. Gouesbet, B. Maheu, and G. Grehan, Journal of the Optical Society of America A 5, 1427 (1988).
  * Huang (2013) K. Huang, Master’s thesis, ShangHai University (2013).
  * Wang (2016) H. Wang, Master’s thesis, ShangHai University (2016).
  * Plimpton (1995) S. Plimpton, Journal of Computational Physics 117, 1 (1995).
  * Masayoshi (2009) M. Masayoshi, _Technical Foundation of Film Preparation (4th edition, in Japanese)_ (Chemical Industry Press, 2009).
  * Noguchi et al. (2007) H. Noguchi, N. Kikuchi, and G. Gompper, Epl 78 (2007).
  * Groot and Warren (1997) R. D. Groot and P. B. Warren, Journal of Chemical Physics 107, 4423 (1997).
  * Israelachvili (1985) J. N. Israelachvili, _Intermolecular and Surface Forces (Third Edition)_ (Academic Pr.*, 1985).
  * Hamaker (1937) H. C. Hamaker, Physica 4, 1058 (1937).
  * Goldhirsch (2003) I. Goldhirsch, Annual Review of Fluid Mechanics 35, 267 (2003).
  * Silbert et al. (2002) L. E. Silbert, D. Ertaş, G. S. Grest, T. C. Halsey, and D. Levine, Physical Review E Statistical Nonlinear & Soft Matter Physics 65, 031304 (2002).
  * Brilliantov et al. (1996) N. V. Brilliantov, F. Spahn, J. M. Hertzsch, and T. Poschel, Physical Review E 53, 5382 (1996).
  * Frenkel and Smit (1997) D. Frenkel and B. Smit, Physics Today 50, 66 (1997).
  * Ohring (2001) M. Ohring, _Materials Science of Thin Films (Second Edition)_ (2001).
  * Zhao and Luo (2005) K. H. Zhao and W. Y. Luo, _Thermology_ (Higher Education Press, 2005).
  * Lun (2006) C. Lun, Journal of Fluid Mechanics 233, 539 (2006).
  * Drazin and Reid (1984) P. G. Drazin and W. H. Reid, Physics Today 37, 87 (1984).
  * Moler and Stewart (1973) C. B. Moler and G. W. Stewart, Siam Journal on Numerical Analysis 10, 241 (1973).
  * Nathanson (2004) G. M. Nathanson, Annual Review of Physical Chemistry 55, 231 (2004).
  * Anderson and Fenn (2004) J. B. Anderson and J. B. Fenn, Physics of Fluids 8, 780 (2004).
  * Moseler and Landman (2000) M. Moseler and U. Landman, Science 289, 1165 (2000).

###### List of Figures

  1. 1  Conceptual prototype of laser induced deposition. . Simplified computational model as particles emitting from a point source into a Gauss laser beam.
  2. 2  Comparison of the distribution and theoretical distribution of the crossing density at the counting plane $h=10$ under different power influences. . Density distribution at different powers. . Non-dimensionalized distribution.
  3. 3 Hydrodynamic stability analysis of light forces on particle jets. Normalized basic state radial force $\bar{F}_{r}$ and particle temperature $\bar{T}$. Effect of electro-Euler number $Eu$ on flow stability.
  4. 4 screenshots of the spatial distribution of particle jets at different powers. (a). $P=0.0W$. (b). $P=0.5W$. (c). $P=1.0W$.
  5. 5 Trajectories of $50$ adjacently numbered particles at $P=1.0W$.
  6. 6 Non-dimensional cumulative distributions of particle jeta at $h=80\mu m$ counting plane at different powers.
  7. 7 Dispersive distribution variations on $h=80\mu m$ counting plane according to changes of electro-Euler number.

###### Contents

  1. I Introduction
  2. II Physical model and simulation methods
    1. II.1 Physical model
    2. II.2 Rayleigh particles affected by laser force
    3. II.3 Numerical aspects
  3. III Different particle source simulations and results
    1. III.1 Point source simulation
    2. III.2 Theoretical analysis and simulation of light converging nanoparticle jets
      1. III.2.1 Light Convergence Particle Jet Stability
      2. III.2.2 Particle Dynamics Simulation of Nanoparticle Jet
  4. IV Conclusion