# Laser Powered Harvesting System for Table-Top Grown Strawberries

Mohamed Sorour, Pål Johan From Robotics Group, Faculty of Science and
Technology, Norwegian University of Life Sciences (NMBU), 1432 Ås, Norway.
<EMAIL_ADDRESS>

###### Abstract

In this paper, a novel tool prototype for harvesting table-top grown
strawberries is presented. With robustness against strawberry localization
error of $\pm 15mm$ and average cycle time of $8.02$ seconds at $50$% of
maximum operational velocity, it provides a promising contribution towards
full automation of strawberry harvesting. In addition, the tool has an overall
fruit-interacting width of $35mm$ that greatly enhances reach-ability due to
the minimal footprint. A complete harvesting system is also proposed that can
be readily mounted to a mobile platform for field tests. An experimental
demonstration is performed to showcase the new methodology and derive relevant
metrics.

###### Index Terms:

fruit picking, strawberry harvesting, agriculture robotics.

##  I Introduction

On account of the seasonal nature of manual harvesting [1, 2], increasing
demand for labours [3, 4] in competitive industrial sectors, as well as ageing
farming society [5], the need for harvesting automation is well established
and justified. Such unmet need, results in elevated cost of living [6] due to
the higher costs incurred to secure labour, and lesser supply of crop, damaged
for missing the optimal harvesting window. Harvesting robots developed so far,
on the other hand, have an average success rate as low as $66\%$ [7], subject
to drastic drop in cluttered environment, mostly due to bulky harvesting tools
[8].

Figure 1: Harvesting setup simulating an actual table-top strawberry growing
system in poly-tunnels.

Automated fruit harvesting essentially requires a tool that can: (1) capture
the fruit, and (2) detach it. To capture, grasping is conventionally employed
in the literature, harvesting apples [9, 10, 11], plums [12], kiwi [13],
tomatoes [14], sweet pepper [15, 16], and strawberries [17]. Mechanical
grasping mostly results in a bulky tool interacting with the produce,
hindering as such reach-ability. Employing vacuum suction instead [18, 19, 20,
21, 22], reduces the tool size, but lead to complexity in finding the correct
spot on the produce to apply, as well as being force-full. To detach, sharp
mechanical cutter is by far the most commonly used method, enlarging the
harvesting tool size, since the driving source (motor) and transmission must
be very close to the cutting tool. Non-conventional cutting means are rare in
the literature, featuring an oscillating blade to cut sweet pepper in [23],
and a thermal cutting device in [24, 25] for cucumber and pepper respectively.
Laser beam is used in [26, 27] to showcase the potential for cutting tomato
peduncles, and in [28, 29, 30, 31] for weed control. Despite the non-
conventional cutting, bulky grasping is employed in the aforementioned
research work.

Figure 2: The harvesting system anatomy in side and top views. Relevant
coordinate frames are shown with x, y, and z axes in red, green and blue
arrows respectively. Vision cone of first and second RGB-D cameras in grey and
light red respectively. Selected dimensions of interest in centimeters. Figure
3: The harvesting tool anatomy with the process observing camera omitted for
convenience.

In this work, the authors present a novel harvesting tool that virtually
captures the fruit by surrounding it, and detaches the stem by applying a
highly focused laser beam from large distance, resulting in minimal
interaction with the fruit and its local environment. In addition, the authors
propose the harvesting system prototype shown in Fig. 1, with relative
distances provided in Fig. 2 replicating those normally found in strawberry
growing poly-tunnel setups. The system as such, can be readily mounted to a
mobile robot for field harvesting, with minimal modifications required. The
contribution of our approach is threefold:

  * •

Productive: with average cycle time of $8$ seconds at $50$% of maximum robot
velocity, in addition to the maintenance-free laser cutting.

  * •

Small footprint: interacting hardware width of $35mm$ greatly enhances fruit
reachability.

  * •

Robust: precise stem entrapment tolerating strawberry localization error of up
to $\pm 15mm$.

This paper is organised as follows, section II introduces the anatomy of the
harvesting system and the operation logic. Algorithms for strawberry
localization and harvesting as well as experiments are reported in section
III. Conclusions are finally given in section IV.

##  II System Anatomy and Operation Logic

The developed harvesting tool is shown in Fig. 3, two components interact with
the strawberry namely; the stem trapping groove, and the stem trapper. Both
are made of mild steel due to its low thermal conductivity (as compared to
Aluminium for example), however, any ferrous metal will be as effective. When
the laser beam is activated during the stem cutting process, this feature
preserves the heat locally, helping the cutting process, killing plant
viruses, and preventing the melt down of other plastic components of the tool.
The width of both components is $35$ and $30mm$ respectively, matching the
dimensions of a large strawberry, as such, minimum fruit dislocation can be
achieved during harvesting, as well as better fruit reach-ability. The opening
between both parts acts as the strawberry entry side in [32], which will
approach the fruit from below, completely surround it, before performing the
stem cut. A convex lens with focal length of $25cm$ focuses the fiber laser
beam generated at the laser head below [33]. It is servo controlled to perform
a lateral, reciprocating, straight line motion to drive the focal point into a
line with controllable length to effectively cut the trapped strawberry-stem.
Two infrared photo interrupters are used to detect the free falling strawberry
after being detached, that is then led to a storage unit by means of a smooth
plastic passage, protecting the fruit surface from damage.

The harvesting tool is mounted to a cost-effective $6$ DOF collaborative robot
arm system as shown in Fig. 2 in side and top views. A set of two RGB-D
cameras mounted in different view points facing the strawberry growing trough,
to the front and from below. The objective of such arrangement is to obtain a
realistic point cloud of the scene with enough information to determine the
location and the size of strawberries. The field of view of both cameras is
depicted in the same figure, colored in light grey and red respectively,
showing the common area of focus, at which strawberries can be mutually
detected. This setup is based on an an average table-top height of $103cm$
readily available at the poly-tunnel of the host institution, and the average
Thorvald mobile robot height [34] widely used in field robotics. The robot arm
base frame $\bm{\mathcal{F}}_{b}$ is related to the first
$\bm{\mathcal{F}}_{c_{1}}$ and second $\bm{\mathcal{F}}_{c_{2}}$ camera frames
by the constant transformation matrices ${{}^{b}\bm{\textbf{T}}_{c_{1}}}$, and
${{}^{b}\bm{\textbf{T}}_{c_{2}}}$ respectively, and rigidly connected via the
arm base plate. As such, the setup (arm, tool, and cameras) can be fitted to a
mobile robot for field testing with minimal modification to the control
algorithm. The tool tip frame $\bm{\mathcal{F}}_{t}$ is controlled to surround
the fruit prior to harvesting, by manipulating the arm, it aligns with the tip
of the stem trapping groove shown in Fig. 3. Let the pose
${}^{b}\textbf{p}_{t}=\begin{bmatrix}{}^{b}\textbf{r}_{t}^{\top}&\\!\\!\\!{}^{b}\bm{\theta}_{t}^{\top}\end{bmatrix}^{\top}\\!\\!\\!\in\mathbb{R}^{6}$
of the manipulator’s tool tip frame $\bm{\mathcal{F}}_{t}$ expressed in the
base frame $\bm{\mathcal{F}}_{b}$ define the task space coordinates. With
${}^{b}\textbf{r}_{t}=\begin{bmatrix}x&y&z\end{bmatrix}^{\top}$ being the
position vector, and
${}^{b}\bm{\theta}_{t}=\begin{bmatrix}\gamma&\beta&\alpha\end{bmatrix}^{\top}$
denoting a minimal representation of orientation (roll, pitch, yaw RPY
variation of Euler angles). In this work, the tool tip frame orientation is
fixed and identical to that of the base frame, as such, the tool pose
${}^{b}\textbf{p}_{t}$ notation will be dropped, and only the $3D$ tool
position ${}^{b}\textbf{r}_{t}$ will be used in the sequel. The arm
configuration shown in Fig. 2 corresponds to the HOME tool position
${}^{b}\textbf{r}_{\texttt{HOME}}$, at which, the arm is not obstructing the
field of vision of the cameras, and the strawberry localization algorithm (to
be presented in the next section) can run.

The operation logic of the harvesting tool is depicted in Fig. 4. Following
strawberry localization, the v-shaped stem-trapping-groove (note the sliding
trapper is also v-shaped, but inverted) - refer to Fig. 4(c)(left) - is
positioned immediately behind the berry-to-harvest at a lower z-axis position
in Fig4(a), the tool tip then elevates to isolate the fruit from its
surrounding to the rear as shown in Fig. 4(b). The trapper then slides to trap
the stem, with Fig. 4(c)(left) showing the first contact between both, and the
eventually fully trapped stem in 4(c)(right). The laser beam is then activated
to full power as in Fig. 4(d), and the convex lens reciprocates laterally with
a stroke of $6mm$ driving the focal point across the stem back and forth until
the stem is cut. The mid-stroke is adjusted to the groove center point.
Although the current tool design can force the stem into a precise location as
compared to the author’s previous work [32], a laser beam with large focal
distance is still needed to minimize the volume of the hardware interacting
with the fruit during the harvesting process. The trapper end of stroke is
adjusted to minimize the stem relocation, and as such the strawberry, while
maintaining a metallic background shield as the laser beam performs the cut,
this is shown in the enlarged square in 4(d).

Figure 4: Operation logic, with the harvesting tool at the x-y coordinates of
the strawberry in (a), eventually encapsulating it by moving upwards in z-axis
(b). The stem is precisely entrapment into a groove in (c) while tolerating
fruit localization error, followed by triggering the laser beam in (d) until
fruit detachment.

The tool, aims to approach the berry from below while perfectly aligned with
the trapper groove. However, due to inaccuracies resulting from the point
cloud calculation, and augmenting the depth information from two different
sensors, localization errors can add up, this in addition to the naturally
occurring stem bending. The width of the trapper of $30mm$ ensures robustness
against such errors within a tolerance of $\pm{15mm}$ from the actual
strawberry location. This is shown in Fig. 4(b) and Fig. 4(c)(left), with
imperfect strawberry localization. Fruit detachment is detected using two
photo interrupters, the infrared beam of which is virtually shown in red color
in Fig. 4(b) and 4(c) for clarification. Once any of the two beams is
interrupted, the laser source is disengaged, terminating the cutting cycle,
and the arm either moves to the next fruit to cut, or to the home position.

Figure 5: Reduced scene point cloud focusing on the area of interest with
black bounding box marking the located strawberries. Algorithm 1 Strawberry
localization algorithm

0: colored scene point cloud ${}^{c_{1}}\mathcal{C}_{s_{1}}$ from RGB-D CAM1,
colored scene point cloud ${}^{c_{2}}\mathcal{C}_{s_{2}}$ from RGB-D CAM2.

0: Vectors defining the bounding box for localized strawberries in arm base
frame.

1:
${}^{b}\mathcal{C}_{s}={{}^{b}\bm{\textbf{T}}_{c_{1}}}^{c_{1}}\mathcal{C}_{s_{1}}+{{}^{b}\bm{\textbf{T}}_{c_{2}}}^{c_{2}}\mathcal{C}_{s_{2}}$

2: ${}^{b}\mathcal{C}_{rs}=\emptyset$

3: for each point ${}^{b}{c}^{i}_{s}$ in ${}^{b}\mathcal{C}_{s}$ do

4: if (${}^{b}{c}^{i}_{s}.x<x^{+}_{l}$ and ${}^{b}{c}^{i}_{s}.x>x^{-}_{l}$ and
${}^{b}{c}^{i}_{s}.y<y^{+}_{l}$ and ${}^{b}{c}^{i}_{s}.y>y^{-}_{l}$ and
${}^{b}{c}^{i}_{s}.z<z^{+}_{l}$ and ${}^{b}{c}^{i}_{s}.z>z^{-}_{l}$ ) then

5: ${}^{b}\mathcal{C}_{rs}\leftarrow{{}^{b}{c}^{i}}_{s}$

6: end if

7: end for

8: ${}^{b}\mathcal{C}_{red}=\emptyset$

9: for each point ${}^{b}{c}^{i}_{rs}$ in ${}^{b}\mathcal{C}_{rs}$ do

10: if (${}^{b}{c}^{i}_{rs}.r>r_{th}$ and ${}^{b}{c}^{i}_{rs}.g<g_{th}$ and
${}^{b}{c}^{i}_{rs}.b<b_{th}$ ) then

11: ${}^{b}\mathcal{C}_{red}\leftarrow{{}^{b}{c}^{i}_{rs}}$

12: end if

13: end for

14:
${}^{b}\mathcal{C}_{straw}=\texttt{EuCS(}t,s_{min},s_{max},^{b}\mathcal{C}_{red}\texttt{)}$

15: for each cluster ${}^{b}\textbf{c}^{i}_{straw}$ in
${}^{b}\mathcal{C}_{straw}$ do

16:
${}^{b}\textbf{x}_{straw}^{min}\leftarrow\texttt{MIN(}^{b}\textbf{c}^{i}_{straw},\texttt{0)}$

17:
${}^{b}\textbf{x}_{straw}^{max}\leftarrow\texttt{MAX(}^{b}\textbf{c}^{i}_{straw},\texttt{0)}$

18:
${}^{b}\textbf{y}_{straw}^{min}\leftarrow\texttt{MIN(}^{b}\textbf{c}^{i}_{straw},\texttt{1)}$

19:
${}^{b}\textbf{y}_{straw}^{max}\leftarrow\texttt{MAX(}^{b}\textbf{c}^{i}_{straw},\texttt{1)}$

20:
${}^{b}\textbf{z}_{straw}^{min}\leftarrow\texttt{MIN(}^{b}\textbf{c}^{i}_{straw},\texttt{2)}$

21:
${}^{b}\textbf{z}_{straw}^{max}\leftarrow\texttt{MAX(}^{b}\textbf{c}^{i}_{straw},\texttt{2)}$

22: end for

23: return
${}^{b}\textbf{x}_{straw}^{min},{{}^{b}\textbf{x}_{straw}^{max}},{{}^{b}\textbf{y}_{straw}^{min}},{{}^{b}\textbf{y}_{straw}^{max}},$
${}^{b}\textbf{z}_{straw}^{min},{{}^{b}\textbf{z}_{straw}^{max}}$

Figure 6: Snapshots of a complete strawberry harvesting cycle, starting and
ending with the detachment of a single strawberry, indicating the time stamp
(in seconds) of each step in-between. Lower row showing the close-up images
captured by the tool-observing camera.

##  III Control Algorithm and Experiments

In this section we present the algorithms for strawberry localization and
motion control that are used in the harvesting experiment. In which, the
system described in the previous section is used to harvest a collection of
$9$ strawberries in full autonomy, the point cloud of these depicted in Fig.
5. The hardware used feature the $6$ DOF xARM6 collaborative robot [35], two
realsense D435 depth cameras, and a $50$ Watts Raycus RFL-P50QB module [33] as
the fiber laser source. The harvesting tool is 3D printed except for the stem-
trapper and the trapping-groove that are subjected to extensive laser heat.
Microcontroller is used to control the convex lens and trapper movements, the
laser activation, as well as monitoring the strawberry detachment photo
interrupters. On the software side, we use the realsense SDK library for
interfacing with the cameras, and the Point Cloud Library (PCL) [36]. A video
of the harvesting experiment is submitted with this paper, also available on
YouTube111https://youtu.be/W3UyDt_7erA., snapshots of which are shown in Fig.
6.

###  III-A Strawberry Localization

The pseudo-code for strawberry localization is provided in Algorithm 1, it
takes as input the raw scene point cloud acquired by both RGB-D cameras
${}^{c_{1}}\mathcal{C}_{s_{1}}$, and ${}^{c_{2}}\mathcal{C}_{s_{2}}$ expressed
in the corresponding camera frame. These are then transformed to the arm base
frame and augmented (code line 1) to form the scene point cloud
${}^{b}\mathcal{C}_{s}$. From which, a reduced scene cloud set
${}^{b}\mathcal{C}_{rs}\subset{{}^{b}\mathcal{C}}_{s}$ is then constructed,
where a point ${{}^{b}{c}^{i}}_{s}$ in the scene cloud ${}^{b}\mathcal{C}_{s}$
is added to the reduced cloud if it resides in a spatial window characterized
by the maximum and minimum limits
$x^{+}_{l},x^{-}_{l},y^{+}_{l},y^{-}_{l},z^{+}_{l},z^{-}_{l}$ in $x,y,$ and
$z$ coordinates respectively of the arm base frame. This reduces the
forthcoming computation to the area of interest that is dexterously accessible
by the robot arm, such reduced point cloud is shown in Fig. 5. Ripe
strawberries are simply extracted by thresholding the red color using the RGB
thresholds $r_{th},g_{th},b_{th}$ to form the red point cloud
${}^{b}\mathcal{C}_{red}\subset{{}^{b}\mathcal{C}}_{rs}$. We then use the PCL
implemented Euclidean cluster segmentation algorithm to segment each
individual strawberry. The output of the function
$\texttt{EuCS(}t,s_{min},s_{max},^{b}\mathcal{C}_{red}\texttt{)}$ is a set of
strawberry clusters (set of point clouds each representing a single
strawberry) ${}^{b}\mathcal{C}_{straw}$ arranged in ascending order of the
y-axis coordinate values. With $t,s_{min},s_{max}$ denoting the segmentation
tolerance, minimum and maximum cluster sizes respectively. The
$\texttt{MIN(}^{b}\textbf{c}^{i}_{straw},\texttt{0)}$ function, supplied with
a point cloud ${}^{b}\textbf{c}^{i}_{straw}$ and an index, will return the
minimum value available at such index in all point cloud points, whereas,
supplied with a vector, will return the smallest value irrespective of the
index. The lengthy pseudo-code for the aforementioned two functions is omitted
for convenience. Values of the parameters used in Algorithm 1 is provided in
Table I. The output is a set of vectors defining the bounding boxes of the
localized strawberries, these boxes are shown in Fig. 5.

TABLE I: Localization algorithm parameters Parameter | Value(unit) | Parameter | Value(unit)
---|---|---|---
$x^{+}_{l}$ | $55(cm)$ | $r_{th}$ | $100$
$x^{-}_{l}$ | $25(cm)$ | $g_{th}$ | $70$
$y^{+}_{l}$ | $30(cm)$ | $b_{th}$ | $70$
$y^{-}_{l}$ | $-30(cm)$ | $t$ | $0.02(cm)$
$z^{+}_{l}$ | $50(cm)$ | $s_{min}$ | $20$
$z^{-}_{l}$ | $30(cm)$ | $s_{max}$ | $1000$

###  III-B Harvesting Motion

Algorithm 2 Harvesting algorithm

0:
${}^{b}\textbf{x}_{straw}^{max},{{}^{b}\textbf{y}_{straw}^{min}},{{}^{b}\textbf{y}_{straw}^{max}},{{}^{b}\textbf{z}_{straw}^{min}},{{}^{b}\textbf{z}_{straw}^{max}}$.

0: desired tool tip position in arm base frame ${}^{b}\textbf{r}_{t}^{d}$.

1: $\texttt{ROBOT\\_MOVE(}^{b}\textbf{r}_{\texttt{HOME}}\texttt{)}$

2:
${}^{b}{z}_{min}=\texttt{MIN(}^{b}\textbf{z}_{straw}^{min}\texttt{)}-\texttt{10mm}$

3: for each strawberry unit i in $n_{straw}$ do

4: ${}^{b}\textbf{r}_{t}^{d}\texttt{(2)}={{}^{b}{z}_{min}}$

5: $\texttt{ROBOT\\_MOVE(}^{b}\textbf{r}_{t}^{d}\texttt{)}$

6:
${}^{b}\textbf{r}_{t}^{d}\texttt{(0)}={{}^{b}\textbf{x}_{straw}^{max}\texttt{(i)}}+\texttt{10mm}$

7:
${}^{b}\textbf{r}_{t}^{d}\texttt{(1)}=({{}^{b}\textbf{y}_{straw}^{min}\texttt{(i)}}+{{}^{b}\textbf{y}_{straw}^{max}\texttt{(i)}})/2$

8: $\texttt{ROBOT\\_MOVE(}^{b}\textbf{r}_{t}^{d}\texttt{)}$

9:
${}^{b}\textbf{r}_{t}^{d}\texttt{(2)}={{}^{b}\textbf{z}_{straw}^{max}\texttt{(i)}}+\texttt{15mm}$

10: $\texttt{ROBOT\\_MOVE(}^{b}\textbf{r}_{t}^{d}\texttt{)}$

11: TRAP_STEM()

12: while IR1 and IR2 do

13: LASER(ON)

14: end while

15: LASER(OFF)

16: RELEASE_STEM()

17: end for

18: $\texttt{ROBOT\\_MOVE(}^{b}\textbf{r}_{\texttt{HOME}}\texttt{)}$

Pseudo-code of the harvesting algorithm is provided in Algorithm 2, it is also
depicted in Fig. 6. It takes as input, the vectors defining the bounding box
of the identified/localized strawberries, that is the output of the
localization algorithm in Algorithm 1. The output is a series of commands to
control the robot arm movement, as well as controlling the stem-trapper and
the laser module. In this work, a strawberry cutting cycle, and as such the
cycle time, starts and ends with a single strawberry unit being detached.
Initially, the algorithm moves the robot arm to the home position
$\textbf{r}_{\texttt{HOME}}$ (configuration depicted in Fig. 2), the
ROBOT_MOVE() function in Algorithm 2 is blocking (has to be finished before
executing the next line of code). It then computes a minimal z-axis point
${}^{b}{z}_{min}$, corresponding to the lowest hanging strawberry observed. In
Fig. 6(a), a strawberry unit has just been cut, signaling the start of a new
cycle, at which the robot arm moves to ${}^{b}{z}_{min}$, as shown in Fig.
6(b) consuming $1.77$ seconds. By setting
${}^{b}\textbf{r}_{t}^{d}\texttt{(2)}={{}^{b}{z}_{min}}$, only the z-axis
component of the desired tool position vector ${}^{b}\textbf{r}_{t}^{d}$ is
changed, while retaining the old values for the other two components (x and y
axes). The arm then moves to the x-y coordinates of the next detected
strawberry using the information of the corresponding bounding box, before
moving upwards in z-axis as shown in Fig. 6(c) and Fig. 6(d) respectively. In
order to compensate for the probable inaccuracy in measuring the actual depth
of strawberry, an empirically determined safety factor is added to the maximum
boundary computed in x and z axes (code lines 6 and 9 in Algorithm 2). To this
end, the robot arm motion for this cycle is finished, what remains is the
cutting procedure, that is microcontroller implemented. When given the signal
to perform the cut, the microcontroller actuates the trapper forward using
TRAP_STEM(), then activates the laser module as well as the convex lens
lateral movement using LASER(ON) until any of the photo interrupters IR1 or
IR2 (refer to Fig. 4) detects the detached fruit. At which instance, the laser
beam is deactivated and the trapper moves backwards.

The average cycle and stem-cut time is $8.02$, and $2.3$ seconds respectively,
with the robot arm operating at $50$% of maximum velocity for hardware safety
reasons. The localization algorithm consumes $100ms$ in worst case scenario on
a standard laptop with intel core-i$7$ processor. The maintenance-free laser
cutting has a long term positive impact on productivity, as compared to
literature-dominating mechanical cutting tools that require frequent
replacement. To this end, the authors believe that a $3D$ linear system
replacing the robot arm, in addition to a $100$ Watt laser module (following
the conclusions in [37]) can reduce the average cycle time below the $4$
seconds mark, cutting the gap with the average manual picking cycle standing
at $1.5$s [38]. That, in addition to further enhancing the fruit
identification and localization algorithm, would greatly update the current
development status [39] thanks to the minimal footprint of the developed tool.

##  IV Conclusion

In this work, a novel harvesting prototype has been presented, customised for
table-top grown strawberries. It employs a trapping mechanism to surround the
fruit and force the stem into a pre-defined position where a focused laser
beam performs the cut. The proposed harvesting tool is robust against
localization errors while maintaining small footprint for enhancing fruit
reachability. Successful harvesting demonstration confirm the effectiveness of
the novel methodology.

## References

  * [1] K. Nolte and M. Ostermeier, “Labour market effects of large-scale agricultural investment: Conceptual considerations and estimated employment effects,” _World Development_ , vol. 98, pp. 430–446, 2017.
  * [2] U. Government, “The impact on the horticulture and food processing sectors of closing the seasonal agricultural workers scheme and the sectors based scheme,” https://assets.publishing.service.gov.uk/government/uploads/system/uploads/attachment_data/file/257242/migrant-seasonal-workers.pdf, last accessed 18/02/2022, May 2013.
  * [3] T. W. Bank, “Employment in agriculture (percentage of total employment) (modeled ilo estimate),” https://data.worldbank.org/indicator/SL.AGR.EMPL.ZS, accessed 12/06/2021.
  * [4] N. F. U. UK, “Establishing the labour availability issues of the uk food and drink sector,” https://www.nfuonline.com/archive?treeid=152097, last accessed 18/02/2022, August 2021.
  * [5] T. Duckett, S. Pearson, S. Blackmore, and B. Grieve, “Agricultural robotics: The future of robotic agriculture,” _EPSRC UK Robotics and Autonomous Systems (RAS) Network_ , 2018. [Online]. Available: http://arxiv.org/abs/1806.06762
  * [6] A. Cassey, K. Lee, J. Sage, and P. Tozer, “Assessing post-harvest labor shortages, wages, and welfare,” _Agricultural and Food Economics_ , vol. 6, 2018.
  * [7] C. W. Bac, E. J. van Henten, J. Hemming, and Y. Edan, “Harvesting robots for high-value crops: State-of-the-art review and challenges ahead,” _Journal of Field Robotics_ , vol. 31, no. 6, pp. 888–911, 2014.
  * [8] G. Kootstra, X. Wang, P. M. Blok, J. Hemming, and E. van Henten, “Selective harvesting robotics: Current research, trends, and future directions,” _Current Robotics Reports_ , vol. 2, pp. 95–104, 2021.
  * [9] Z. De-An, L. Jidong, J. Wei, Z. Ying, and C. Yu, “Design and control of an apple harvesting robot,” _Biosystems Engineering_ , vol. 110, no. 2, pp. 112–122, 2011. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S1537511011001206
  * [10] A. Silwal, J. R. Davidson, M. Karkee, C. Mo, Q. Zhang, and K. Lewis, “Design, integration, and field evaluation of a robotic apple harvester,” _Journal of Field Robotics_ , vol. 34, no. 6, pp. 1140–1159, 2017. [Online]. Available: https://onlinelibrary.wiley.com/doi/abs/10.1002/rob.21715
  * [11] Y. Onishi, T. Yoshida, H. Kurita, T. Fukao, H. Arihara, and A. Iwai, “An automated fruit harvesting robot by using deep learning,” _ROBOMECH Journal_ , vol. 6, 2019.
  * [12] J. Brown and S. Sukkarieh, “Design and evaluation of a modular robotic plum harvesting system utilizing soft components,” _Journal of Field Robotics_ , vol. 38, no. 2, pp. 289–306, 2021.
  * [13] L. Mu, G. Cui, Y. Liu, Y. Cui, L. Fu, and Y. Gejima, “Design and simulation of an integrated end-effector for picking kiwifruit by robot,” _Information Processing in Agriculture_ , vol. 7, no. 1, pp. 58–71, 2020\. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S2214317318304372
  * [14] Q. Feng, X. Wang, G. Wang, and Z. Li, “Design and test of tomatoes harvesting robot,” in _2015 IEEE International Conference on Information and Automation_ , 2015, pp. 949–952.
  * [15] B. Arad, J. Balendonck, R. Barth, O. Ben-Shahar, Y. Edan, T. Hellström, J. Hemming, P. Kurtser, O. Ringdahl, T. Tielen, and B. van Tuijl, “Development of a sweet pepper harvesting robot,” _Journal of Field Robotics_ , vol. 37, no. 6, pp. 1027–1039, 2020.
  * [16] C. W. Bac, J. Hemming, B. van Tuijl, R. Barth, E. Wais, and E. J. van Henten, “Performance evaluation of a harvesting robot for sweet pepper,” _Journal of Field Robotics_ , vol. 34, no. 6, pp. 1123–1139, 2017.
  * [17] Y. Xiong, C. Peng, L. Grimstad, P. J. From, and V. Isler, “Development and field evaluation of a strawberry harvesting robot with a cable-driven gripper,” _Computers and Electronics in Agriculture_ , vol. 157, pp. 392–402, 2019. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0168169918312456
  * [18] J. Baeten, K. Donné, S. Boedrij, W. Beckers, and E. Claesen, _Autonomous Fruit Picking Machine: A Robotic Apple Harvester_. Berlin, Heidelberg: Springer Berlin Heidelberg, 2008, pp. 531–539.
  * [19] K. Tanigaki, T. Fujiura, A. Akase, and J. Imagawa, “Cherry-harvesting robot,” _Computers and Electronics in Agriculture_ , vol. 63, no. 1, pp. 65–72, 2008, special issue on bio-robotics. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0168169908000458
  * [20] S. Hayashi, K. Shigematsu, S. Yamamoto, K. Kobayashi, Y. Kohno, J. Kamata, and M. Kurita, “Evaluation of a strawberry-harvesting robot in a field test,” _Biosystems Engineering_ , vol. 105, no. 2, pp. 160–171, 2010. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S1537511009002797
  * [21] S. Hayashi, S. Yamamoto, S. Tsubota, Y. Ochiai, K. Kobayashi, J. Kamata, M. Kurita, H. Inazumi, and R. Peter, “Automation technologies for strawberry harvesting and packing operations in japan 1,” _Journal of Berry Research_ , vol. 4, 2014.
  * [22] G. Hu, C. Chen, J. Chen, L. Sun, A. Sugirbay, Y. Chen, H. Jin, S. Zhang, and L. Bu, “Simplified 4-dof manipulator for rapid robotic apple harvesting,” _Computers and Electronics in Agriculture_ , vol. 199, p. 107177, 2022.
  * [23] C. Lehnert, A. English, C. McCool, A. W. Tow, and T. Perez, “Autonomous sweet pepper harvesting for protected cropping systems,” _IEEE Robotics and Automation Letters_ , vol. 2, no. 2, pp. 872–879, 2017.
  * [24] E. van Henten, J. Hemming, B. van Tuijl, B. van Tuijl, J. Kornet, J. Meuleman, J. Bontsema, and E. van Os, “An autonomous robot for harvesting cucumbers in greenhouses,” _Autonomous Robots_ , vol. 13, no. 3, pp. 241–258, 2002.
  * [25] S. BACHCHE and K. OKA, “Performance testing of thermal cutting systems for sweet pepper harvesting robot in greenhouse horticulture,” _Journal of System Design and Dynamics_ , vol. 7, no. 1, pp. 36–51, 2013.
  * [26] J. Liu, Z. Li, P. Li, and H. Mao, “Design of a laser stem-cutting device for harvesting robot,” in _2008 IEEE International Conference on Automation and Logistics_ , 2008, pp. 2370–2374.
  * [27] J. Liu, Y. Hu, X. Xu, and P. Li, “Feasibility and influencing factors of laser cutting of tomato peduncles for robotic harvesting,” _African Journal of Biotechnology_ , vol. 10, no. 69, pp. 15 552–15 563, 2011.
  * [28] T. Heisel, J. Schou, C. Andreasen, and S. Christensen, “Using laser to measure stem thickness and cut weed stems,” _Weed Research_ , vol. 42, no. 3, pp. 242–248, 2002. [Online]. Available: https://onlinelibrary.wiley.com/doi/abs/10.1046/j.1365-3180.2002.00282.x
  * [29] S. K. Mathiassen, T. Bak, S. Christensen, and P. Kudsk, “The effect of laser treatment as a weed control method,” _Biosystems Engineering_ , vol. 95, no. 4, pp. 497–505, 2006. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S1537511006002984
  * [30] G. Coleman, C. Betters, C. Squires, S. Leon-Saval, and M. Walsh, “Low energy laser treatments control annual ryegrass (lolium rigidum),” _Frontiers in Agronomy_ , vol. 2, p. 35, 2021. [Online]. Available: https://www.frontiersin.org/article/10.3389/fagro.2020.601542
  * [31] M. Nadimi, D. Sun, and J. Paliwal, “Recent applications of novel laser techniques for enhancing agricultural production,” _Laser Physics_ , vol. 31, no. 5, 2021.
  * [32] M. Sorour, P. J. From, K. Elgeneidy, S. Kanarachos, and M. Sallam, “Produce harvesting by laser stem-cutting,” in _2022 IEEE International Conference on Automation Science and Engineering_ , 2022, p. to be published.
  * [33] Raycus, “50w q-switched pulse fiber laser,” https://en.raycuslaser.com/products/50w-q-switched-pulse-fiber-laser.html, last accessed 05/09/2022, 2022.
  * [34] L. Grimstad and P. J. From, “Thorvald ii - a modular and re-configurable agricultural robot,” _IFAC-PapersOnLine_ , vol. 50, no. 1, pp. 4588–4593, 2017, 20th IFAC World Congress. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S2405896317314830
  * [35] xARM, “xarm collaborative robot,” https://www.ufactory.cc/xarm-collaborative-robot, last accessed 05/09/2022, 2022.
  * [36] R. B. Rusu and S. Cousins, “3D is here: Point Cloud Library (PCL),” in _IEEE International Conference on Robotics and Automation (ICRA)_ , Shanghai, China, May 9-13 2011.
  * [37] M. Sorour, P. J. From, K. Elgeneidy, S. Kanarachos, and M. Sallam, “Compact strawberry harvesting tube employing laser cutter,” in _2022 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , 2022, pp. 8956–8962.
  * [38] S. Woo, D. D. Uyeh, J. Kim, Y. Kim, S. Kang, K. C. Kim, S. Y. Lee, Y. Ha, and W. S. Lee, “Analyses of work efficiency of a strawberry-harvesting robot in an automated greenhouse,” _Agronomy_ , vol. 10, no. 11, 2020.
  * [39] H. Zhou, X. Wang, W. Au, H. Kang, and C. Chen, “Intelligent robots for fruit harvesting: recent developments and future challenges,” _Precision Agriculture_ , pp. 1573–1618, 2022.


# Directed Flow of Identified Particles in Au + Au Collisions at
$\mbox{$\sqrt{s_{\mathrm{NN}}}$}=200$ GeV

Aihong Tang (for the STAR Collaboration)  Physics Department, Brookhaven
National Laboratory, PO Box 5000, Upton, NY 11973<EMAIL_ADDRESS>
###### Abstract

STAR’s measurements of directed flow ($v_{1}$) for pions, kaons (charged and
$K_{s}^{0}$), protons and antiprotons for Au + Au collisions at
$\mbox{$\sqrt{s_{\mathrm{NN}}}$}=200$ GeV are presented. Negative $v_{1}(y)$
slope is observed for pions, antiprotons, protons and kaons. The proton
$v_{1}(y)$ slope at midrapidity is found extremely small. Sizable difference
is seen between $v_{1}$ of protons and anti-protons in 5-30% central
collisions. Anti-flow can explain the negative slope, however, it has
difficulties explaining the centrality dependence of the difference between
the $v_{1}(y)$ slopes of protons and anti-protons. The $v_{1}$ excitation
function is presented. Comparisons to model calculations (RQMD, uRQMD, AMPT
and QGSM) are made, and it is found that none of the four models can
successfully describe the data.

##  1 Introduction

In relativistic heavy ion collisions, anisotropic flow describes the
collective motion of particles with respect to the reaction plane, and it is
conventionally characterized [1] by Fourier coefficients

$v_{n}=\langle\cos n(\phi-\psi)\rangle$ (1)

where $\phi$ denotes the azimuthal angle of an outgoing particle, $\psi$ is
the orientation of the reaction plane, and $n$ denotes the harmonic. So far
three such coefficients are found finite at RHIC, namely, directed flow
$v_{1}$, elliptic flow $v_{2}$ and 4th order harmonic flow $v_{4}$. This paper
will focus on the directed flow, the first Fourier coefficient. Directed flow
describes the sideward motion of fragments in ultra-relativistic nuclear
collisions, and it carries early information from the collision [2]. The shape
of directed flow at midrapidity is of special interest because it might reveal
a signature of a possible phase transition from normal nuclear matter to a
Quark-Gluon Plasma (QGP) [3, 4, 5]. It is argued that directed flow, as a
function of rapidity, may exhibit a flatness at midrapidity due to a strong,
tilted expansion of the source. Such expansion gives rise to anti-flow [3] or
a 3rd flow [4] component. The anti-flow (3rd flow component) is perpendicular
to the source surface, and is in the direction that is opposite to the
direction of the bouncing-off motion of nucleons. If the tilted expansion is
strong enough, it can even overcome the bouncing-off motion and results in a
negative $v_{1}(y)$ slope at midrapidity, producing a wiggle-like structure.
In [5], it is emphasized that proton $v_{1}$ can serve as a probe for the
first order phase transition, as calculation shows that proton $v_{1}(y)$ does
not show negative slope for a hadronic Equation of State (EoS) without a QGP
phase transition. Note that although the calculation is done with a first
order phase transition at SPS energies, the direct cause of the negative slope
is the strong, tilted expansion, which is also relevant at top RHIC energies.
A wiggle structure is also seen in RQMD [6], and it is attributed to the
baryon stopping and positive space-momentum correlation. In this picture, no
phase transition is needed, and pions and nucleons flow in opposite
directions. To distinguish between baryon stopping and anti-flow associated
with a phase transition, it is desirable to measure the $v_{1}(y)$ for
identified particles and compare the signs of their slope at midrapidity. The
centrality dependence of proton $v_{1}(y)$ may have implications for a
possible first order phase transition. It is expected that in very peripheral
collisions, protons flow in the same direction as spectators. In mid-central
collisions, if there is a phase transition, the proton $v_{1}(y)$ slope at
midrapidity may change sign and become negative. Eventually the slope
diminishes in central collisions due to the symmetry of collisions.

At low energies, it is shown by the E895 collaboration [7] that K${}^{0}_{s}$
has a negative $v_{1}(y)$ slope around midapidity, while lambda and protons
have positive slope [8]. It is explained by a repulsive kaon-nucleon potential
and attractive lambda-nucleon potential. The NA49 collaboration [9] has
measured the $v_{1}$ for pions and protons, and negative $v_{1}(y)$ slope is
observed by the standard event plane method. The three particle correlation
method $v_{1}\\{3\\}$, which is believed to be less sensitive to non-flow
effects, gives negative slope too, but with larger statistical error. At top
RHIC energies, $v_{1}$ has been studied for charged particles by both the STAR
and the PHOBOS collaborations [10, 11, 12, 13]. It is found that $v_{1}$ in
the forward region follows the limiting fragmentation hypothesis well, and
$v_{1}(\eta)$ depends only on the incident energy, but not on the size of the
colliding system at a given centrality. The system size independence of
$v_{1}$ can be explained by Hydrodynamic calculation with an titled initial
condition [14]. The systematic study of $v_{1}$ for identified particles at
RHIC did not begin until recently because it is more challenging for two
reasons, 1) $v_{1}$ for some identified particles (for example, protons) is
much smaller than that of charged particles, thus is more difficult to measure
, 2) it demands more statistics to measure $v_{1}$ for identified particles
other than poins.

##  2 Data Set

In this analysis we used 62 million events for Au + Au collisions at
$\mbox{$\sqrt{s_{\mathrm{NN}}}$}=200$ GeV, all taken by STAR’s minimum bias
trigger during RHIC run 2007. Charged particles are identified by their energy
loss inside STAR’s Time Projection Chamber (TPC) [15], and
$K_{s}^{0}(\rightarrow\pi^{+}\pi^{-})$ are reconstructed by their charged
daughter tracks inside the TPC. Track quality cuts are the same as used in
[16]. The centrality definition of an event was based on the number of charged
global tracks in the TPC with track quality cuts: $|\eta|<$0.5, a Distance of
Closest Approach (DCA) to the global vertex less than 3 cm, and 10 or more fit
points. Additional weight is assigned to each event in the analysis accounting
for the bias in the vertex Z direction, due to the inclusion of hits of
Silicon Vertex Tracker in tracking for year 2007 data production. In addition,
the transverse momentum $p_{T}$ for protons are required to be larger than 400
MeV$/c$, and DCA are required to be less than 1 cm, in order to avoid the
beam-pipe background. The same cut is applied to anti-protons as well to
enssure a fair comparison with protons. The high end of the $p_{T}$ cut is 1
GeV$/c$ where protons and pions have the same energy loss in the TPC and thus
become indistinguishable. The event plane angle is determined from the
sideward deflection of spectator neutrons measured by STAR’s Shower Maximum
Detector at Zero Degree Calorimeters (ZDCSMD), which are located close to beam
rapidity and have the minimum contribution from non-flow. The description of
measuring $v_{1}$ using the ZDCSMD event plane can be found in [12].

Figure 1: Pion $v_{1}$ as a function of rapidity (or $\eta$, for the QGSM
model only) for 10-70% Au + Au collisions at
$\mbox{$\sqrt{s_{\mathrm{NN}}}$}=200$ GeV. Also plotted are corresponding
model calculations. The calculation from UrQMD and QGSM are scaled by a factor
of 2.5 and 0.18, respectively. Errors are statistical only.

##  3 Result

Fig. 1 shows the pion $v_{1}$ as a function of rapidity. The result is
compared to four model calculations, namely, RQMD [6], UrQMD [17], AMPT [18]
and QGSM [19]. Errors in this plot (and other data figures in this paper) are
statistical only. Systematical uncertainties are studied by varying track
quality cuts, and are found to be negligible. Calculations from RQMD and QGSM
are scaled by a factor of 2.5 and 0.18, respectively. Following convention,
the sign of spectator $v_{1}$ in the forward region is chosen to be positive,
to which the measured sign of $v_{1}$ for particles of interest is only
relative. The pion $v_{1}(y)$ result follows closely to the published charged
particle result [12] and it shows a negative slope at midrapidity. Most models
either predict the wrong sign of pion $v_{1}$ or the wrong magnitude by a
signifiant factor 111However, in a recent paper [14] posted after the
$26^{th}$ Winter Workshop on Nuclear Dynamics, it is shown that the directed
flow of charged particles can be explained by Hydrodynamic calculation with a
tilted initial condition., which indicates that it is a challenging job to
understand the dynamics happening at early stages. In this sense the
measurement presented in this paper will offer strong constraints on models.

Figure 2: Proton (solid circles) and anti-proton (solid squares) $v_{1}$ as a
function of rapidity for 10-70% Au + Au collisions at
$\mbox{$\sqrt{s_{\mathrm{NN}}}$}=200$ GeV, pion $v_{1}(y)$ (empty stars) is
shown as a reference.

To address the shadowing effect and the anti-flow, in Fig.2, proton and anti-
proton $v_{1}$ are plotted together with pion $v_{1}$. Pions are the dominant
particles produced and they carry the flow of the bulk. It is not a surprise
to see that anti-protons, also as produced particles, follow the flow
direction of pions. They both have negative slope at midrapidity, which is
consistent with the anti-flow picture. $v_{1}(y)$ of protons, however,
exhibits a flatter shape than others. The mass difference between protons and
pions cannot explain the flatness of proton $v_{1}(y)$ because the anti-
protons have the same mass as that of protons, yet they flow close to pions.
Indeed, the observed $v_{1}$ for protons is a convolution of directed flow of
produced protons with that of transported protons, and the flatness of
inclusive proton $v_{1}(y)$ around midrapidity can be explained by the
negative flow of produced protons being compensated by the positive flow of
protons transported from spectator rapidity, as a feature expected in the
anti-flow picture.

Figure 3: Charged kaon (up triangles) and $K_{s}^{0}$ (down triangles) $v_{1}$
as a function of rapidity for 10-70% Au + Au collisions at
$\mbox{$\sqrt{s_{\mathrm{NN}}}$}=200$ GeV, pion $v_{1}(y)$ (empty stars) is
shown as a reference.

Fig. 3 shows $v_{1}$ of kaon and $K_{s}^{0}$ as a function of rapidity, for
centrality 10-70%. Pion $v_{1}$ is also plotted for comparison. Kaon $v_{1}$
is a unique measurement because the kaon/proton cross section is smaller than
the pion/proton cross section, thus its directed flow is less likely to be
”induced” by the shadowing effect. The plot shows that both kaon and
$K_{s}^{0}$ have negative $v_{1}$ slope, which means that the observed $v_{1}$
is not likely to be caused by the shadowing effect.

On the other hand, anti-flow has difficulties in explaining the centrality
dependence of $v_{1}$, as shown by Fig. 4 , in which $v_{1}$ slope at
midrapidity is plotted as a function of centrality for protons, anti-protons
and charged particles.

Figure 4: Charged (solid stars), proton (solid circles) and anti-proton (solid
squares) $v_{1}(y)$ slope ($dv_{1}/dy$) at midrapidity as a function of
centrality for Au + Au collisions at $\mbox{$\sqrt{s_{\mathrm{NN}}}$}=200$
GeV.

If there is an anti-flow due to the strong, tilted expansion, one expects such
an effect is larger in mid-central collisions than that in peripheral
collisions. As a consequence, proton $v_{1}$ slope, which is expected to be
positive in very peripheral collisions, will change its sign to negative in
mid-central collisions and approach zero in central collisions. Due to the
large statistical error, we cannot exclude the possibility that proton v1 can
be positive in 70-80% central collisions. However, in 30-80% central
collisions, proton $v_{1}$ slope is found mostly negative and the magnitude
decreases with decreasing centrality. In more central (5-30%) collisions,
proton $v_{1}$ slope becomes extremely small, while anti-proton $v_{1}$ slope
remains negative and continues to follow that of charged particles (mostly
pions). Anti-flow may cause a difference between $v_{1}$ of protons and anti-
protons, and such difference is expected to be accompanied by strongly
negative $v_{1}$. In data, the large difference between proton and anti-proton
$v_{1}$ is seen in 5-30% centrality, while strongly negative $v_{1}$ is found
for protons, anti-protons and charged particles in a different centrality
(30-70%).

Figure 5: Proton $v_{1}(y^{{}^{\prime}})$ slope ($dv_{1}/dy^{{}^{\prime}}$) at
midrapidity as a function of center of mass collision energy, where
$y^{{}^{\prime}}=y/y_{beam}$.

In Fig. 5, proton $v_{1}(y/y_{beam})$ slope at midrapidity is plotted as a
function of collision energy. We see that proton $v_{1}$ slope decreases
rapidly with increasing energy, reaching zero around
$\mbox{$\sqrt{s_{\mathrm{NN}}}$}=9$ GeV. It changes its sign to negative as
shown by the data point at $\mbox{$\sqrt{s_{\mathrm{NN}}}$}=17$ GeV, measured
by NA49 experiment. With previous measurements, which include only one point
above $\mbox{$\sqrt{s_{\mathrm{NN}}}$}=9$ GeV, one cannot conclude if proton
$v_{1}$ slope continues to decrease or stays close to zero when energy
increases, the addition of the data point from RHIC indicates that proton
$v_{1}$ slope remains close to zero at $\mbox{$\sqrt{s_{\mathrm{NN}}}$}=200$
GeV. Judging over the broad energy range, the transition of the trend happens
around $\mbox{$\sqrt{s_{\mathrm{NN}}}$}=9$ GeV, making it a discontinuous
point, and interestingly, it coincides with the energy vicinity where $\langle
k^{+}\rangle/\langle\pi^{+}\rangle$ exhibits a kink [20]. Hopefully the on-
going Beam Energy Scan program will map out the region of interest with
details.

##  4 Summary

For 10-70% central Au + Au collisions at $\mbox{$\sqrt{s_{\mathrm{NN}}}$}=200$
GeV, directed flow for pions, kaons (charged and $K_{s}^{0}$), and anti-
protons are found to have negative slope at mid-rapidity. Protons, however,
exhibit a flat shape around mid-rapidity. Sizable difference is seen between
$v_{1}$ of protons and anti-protons in 5-30% central collisions. Anti-flow can
explain the negative slope, however, it has difficulties explaining the
centrality dependence of the difference between the $v_{1}(y)$ slopes of
protons and anti-protons.

## References

## References

  * [1] A. M. Poskanzer and S. A. Voloshin, Phys. Rev. C  58, 1671 (1998).
  * [2] H. Sorge, Phys. Rev. Lett. 78, 2309 (1997).
  * [3] J. Brachmann, S. Soff, A. Dumitru, H. Stöcker, J. A. Mahruhn, W. Greiner, L. V. Bravina, and D. H. Rischke, Phys. Rev. C 61, 024909 (2000).
  * [4] L. P. Csernai and D. Röhrich, Phys. Lett. B , 458, 454 (1999).
  * [5] H. Stöcker, N͡ucl. Phys. A 750, 121 (2005).
  * [6] R. J. Snellings, H. Sorge, S. A. Voloshin, F. Q. Wang and N. Xu, Phys. Rev. Lett 84, 2803 (2000).
  * [7] P. Chung et al. [e895 Collaboration], Phys. Rev. Lett. 85, 940 (2000).
  * [8] P. Chung et al. [e895 Collaboration], Phys. Rev. Lett. 86, 2533 (2001).
  * [9] C. Alt et al. [na49 Collaboration], Phys. Rev. Lett. 68, 034903 (2003).
  * [10] J. Adams et al. [STAR Collaboration], Phys. Rev. Lett. 92, 062301 (2004).
  * [11] B. B. Back et al. [phobos Collaboration], Phys. Rev. Lett 97, 012301 (2006).
  * [12] J. Adams et al. [STAR Collaboration], Phys. Rev. C  73, 034903 (2006).
  * [13] B. I. Abelev et al. [STAR Collaboration], Phys. Rev. Lett. 101 252301 (2008).
  * [14] Piotr Bożek and Iwona Wyskiel, arXiv:1002.4999 [nucl-th] (2010).
  * [15] M. Anderson et al., Nucl. Instrum. Method. A  499 659 (2003).
  * [16] J. Adams et al. [STAR Collaboration], Phys. Rev. C 72 014904 (2005).
  * [17] M. Bleicher and H. Stöcker, Phys. Rev. B 526 309 (2002).
  * [18] J. Y. Chen, J. X. Zuo, X. Z. Cai, F. Liu, Y. G. Ma and A. H. Tang, Phys. Rev. C 81, 014904 (2010).
  * [19] G. Burau, J. Bleibel, C. Fuchs, A. Faessler, L. V. Bravina, and E. E. Zabrodin, Phys. Rev. C 71 054905 (2005).
  * [20] S. V. Afanasiev et al. [NA49 Collaboration], Phys. Rev. C 66, 054902 (2002); C. Alt et al. [NA49 Collaboration], Phys. Rev. C 77, 024903 (2008); Phys. Rev. C 73, 044910 (2006); T. Anticic et al. [NA49 Collaboration], Phys. Rev. C 69, 024902 (2004).


††thanks: Email<EMAIL_ADDRESS>
# Intense high-order harmonic generation in giant fullerene molecule C240

H.K. Avetissian  Centre of Strong Fields Physics at Research Institute of
Physics, Yerevan State University, Yerevan 0025, Armenia  S. Sukiasyan  Centre
of Strong Fields Physics at Research Institute of Physics, Yerevan State
University, Yerevan 0025, Armenia  Max-Planck-Institut für Kernphysik,
Saupfercheckweg 1, 69117 Heidelberg, Germany  T.M. Markosyan  Centre of Strong
Fields Physics at Research Institute of Physics, Yerevan State University,
Yerevan 0025, Armenia  G.F. Mkrtchian  Centre of Strong Fields Physics at
Research Institute of Physics, Yerevan State University, Yerevan 0025, Armenia

###### Abstract

In this work the extreme nonlinear optical response of a giant fullerene
molecule C240 in strong laser field is studied. The investigation of high-
order harmonic generation in such quantum nanostructure is presented modeling
the C240 molecule and its interaction with the laser field in the scope of the
tight-binding mean-field approach. Electron-electron interaction is modeled by
the parametrized Ohno potentail, which takes into account long-range Coulomb
interaction. The essential role of many body Coulomb interaction in
determining of harmonics intensities is demonstrated. We also consider
vacancy-deffected molecule C240. The presence of a single vacancy breaks the
icosahedral symmetry leading to the emergence of intense even-order harmonics.
We examine the dependence of moderate harmonics on laser frequency that shows
the multiphoton resonant nature of high harmonics generation. The dependence
of cutoff harmonics on both laser intensity and frequency are examined too.

##  I Introduction

Intense light interaction with nanostructures can excite the electrons of the
system through multiphoton channels, leading to extreme nonequilibrium states
[1]. The excited electrons subsequently emit coherent electromagnetic
radiation, encompassing tens to hundreds of harmonics of the incident light
[2, 3]. This fundamental process in intense laser-matter interaction is known
as high harmonic generation (HHG) phenomenon [4, 5]. In atoms, HHG has been
widely used to produce coherent extreme ultraviolet radiation, allowing access
to the extreme time resolution of the underlying quantum processes and
enabling attosecond physics [6, 7]. Among the diverse range of nanostructured
materials suitable for nonlinear extrime optical applications, carbon
allotropes hold a central position [8, 9]. One of the carbon allotropes are
fullerenes [10] which are large molecules formed by closing a graphite sheet,
where the required curvature is achieved by incorporating twelve pentagons
among a given number of graphene hexagons. The most well-known fullerene is
the buckminsterfullerene C60 [11], which possesses icosahedral symmetry. The
discovery of fullerene C60 through laser evaporation of graphite was triggered
the study of many other fullerene molecules. Larger fullerenes, often referred
to as giant fullerenes, can also be constructed with icosahedral symmetry
[12]. These large fullerenes can be visualized as cut-out pieces of graphene
that are folded into an icosahedron. Consequently, they exhibit similar
properties to graphene [13] or graphene quantum dots [14], while remaining
stable due to their closed topological structure. Note that in continuous
limit C60 and related molecules are well described by the Dirac equation in
the curved space and in the field of a monopole [15, 16]. Giant or large
fullerenes have been the subject of active research since the 1990s. For a
more comprehensive overview, we refer the reader to references[17, 18, 19, 20,
21, 22, 23] for earlier studies and references [24, 25, 26, 27, 28, 29, 30]
for more recent investigations.

In the field of HHG, enhancing conversion efficiency is of utmost importance.
This efficiency strongly relies on the density of emitters and the density of
states of these emitters. To this end, molecular systems, clusters, and
crystals have shown potential in significantly increasing harmonic intensity
compared to atomic systems, as they can exploit multiple excitation channels
[31, 32, 33]. As a result, there has been a growing interest in extending HHG
to carbon-based materials, such as semimetallic graphene [34, 35, 36, 37, 38,
39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50], graphene quantum dots [51,
52, 53, 54], and fullerenes [55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66].
. Experimental studies, namely Refs. [59, 60], have reported a robust harmonic
signal from C60 plasma. Additionally, theoretical works have predicted strong
HHG from both C60 [56, 57, 65, 66] and C70 molecules [65] and solid C60 [64].
Notably, the increase in conducting electrons in fullerene molecules leads to
a subsequent rise in density of states, thereby opening up new channels that
can amplify the HHG signal. Consequently, exploring the HHG process in giant
fullerenes becomes a compelling area of interest. With the increasing
fullerene size, the molecules are subject to various types of defects.
Therefore, investigating the impact of defects on HHG in large fullerenes
holds significance. Recent research involve effects of disorder, impurities,
and vacancies on HHG in solids [67, 68, 69, 70, 71, 72, 73, 74, 75]. These
studies have revealed that an imperfect lattice can enhance HHG compared to a
perfect lattice, especially when considering doping-type impurities or
disorders. For C60 and C180, it has been shown that both diagonal and off-
diagonal disorders break inversion symmetry, lift the degeneracy of states,
and create new channels for interband transitions, resulting in enhanced high
harmonic emission [66]. This raises intriguing questions about how vacancies
specifically affect the HHG spectra in large fullerenes. Vacancies can occur
naturally or be introduced in fullerenes through laser or ion/electron
irradiation [76, 77]. Taking into account that vacancy defects introduce
localized electronic states [78] and the HHG process is highly sensitive to
electron wave functions, we can expect new effects in the HHG process at
consideration of vacancy-defected fullerenes.

In this study, we present a microscopic theory that explores the extreme
nonlinear interaction of normal and single vacancy-defected fullerene C240
with strong electromagnetic radiation. Particularly, we consider coherent
interaction with a linearly polarized electromagnetic radiation taking into
account collective electron-electron interactions. Employing the dynamical
Hartree-Fock approximation, we reveal the general and basal structure of the
HHG spectrum and its relation to molecular excitations and icosahedral
symmetry breaking of giant molecules.

The paper is organized as follows. In Sec. II, the model and the basic
equations are formulated. In Sec. III, we present the main results. Finally,
conclusions are given in Sec. IV.

##  II The model and theoretical approach

We start by describing the model and theoretical approach. Fullerene molecule
C240 and C240 with a monovacancy is assumed to interact with a mid-infrared or
visible laser pulse that excites electron coherent dynamics. For the brevity
we refer vacancy-defected C240 molecule as C239. The schematic structure of
these fullerene molecules are deployed in Fig. 1. We assume a neutral
molecules, which will be described in the scope of the tight-binding (TB)
theory. The electron-electron interaction (EEI) is described in the extended
Hubbard approximation [79, 80, 65]. Hence, the total Hamiltonian reads:

$\widehat{H}=\widehat{H}_{0}+\widehat{H}_{\mathrm{int}},$ (1)

where

$\widehat{H}_{0}=-\sum_{\left\langle
i,j\right\rangle\sigma}t_{ij}c_{i\sigma}^{\dagger}c_{j\sigma}+\frac{U}{2}\sum_{i\sigma}n_{i\sigma}n_{i\overline{\sigma}}+\frac{1}{2}\sum_{i,j}V_{ij}n_{i}n_{j}$
(2)

is the free fullerene Hamiltonian. Here $c_{i\sigma}^{\dagger}$ creates an
electron with spin polarization $\sigma=\left\\{\uparrow,\downarrow\right\\}$
at site $i$ ($\overline{\sigma}$ is the opposite to $\sigma$ spin
polarization), and $\left\langle i,j\right\rangle$ runs over all the first
nearest-neighbor hopping sites with the hopping integral $t_{ij}$ between the
nearest-neighbor atoms at positions $\mathbf{r}_{i}$ and $\mathbf{r}_{j}$. The
density operator is: $n_{i\sigma}=c_{i\sigma}^{\dagger}c_{i\sigma}$, and the
total electron density for the site $i$ is:
$n_{i}=n_{i\uparrow}+n_{i\downarrow}$. The second and third terms in Eq. (2)
describe the EEI Hamiltonian, with the parameters $U$ and $V_{ij}$
representing the on-site, and the long-range Coulomb interactions,
respectively. The involved molecules contain single and double carbon bonds,
for which model Hamiltonian (2) has been parameterized extensively over the
years. The input Cartesian coordinates for C240 are obtained from the Yoshida
database [81]. In the present paper, as first approximation, monovacancy is
simulated by removing one carbon atom. The initial structures are further
optimized with the help of IQmol programm [82]. Hence, in the vicinity of the
vacancy the bond lenghts are changed. There is also scenary when the structure
undergoes a bond reconstruction in the vicinity of the vacancy [83]. In either
case, a local distortion of the lattice takes place resulting states that are
strongly localized around defects [84, 85]. For the one-electron hopping
matrix elements, which in this work have been restricted to the nearest
neighbors, we use values close to the graphene hopping matrix elements. The
common choice of hopping matrix element is $t_{0}=2.7$ eV, corresponding to
the C-C bond length of $d_{0}=1.42\mathrm{\mathring{A}}$, while for shorter or
longer bonds, its value is extrapolated using the linear relationship
$t_{ij}=t_{0}+\alpha\left(d_{0}-\left|\mathbf{r}_{i}-\mathbf{r}_{j}\right|\right)$,
with $\alpha=3.5\ \mathrm{eV/\mathring{A}}$ being the electron-phonon coupling
constant. The EEI is modeled by the Ohno potential [86]:

$V_{ij}=\frac{U}{\sqrt{1+\frac{U^{2}\left|\mathbf{r}_{i}-\mathbf{r}_{j}\right|^{2}}{V^{2}d_{m}^{2}}}},$
(3)

where $V$ means the strength of the long range Coulomb interaction, and
$d_{m}$ is the average bond length. Depending on the screening effects a
popular choice of parameters for the Coulomb interactions is $0\leq U\leq
4t_{0}$, and $V=0.5U$ [80, 87].

The light-matter interaction is described in the length-gauge

$\widehat{H}_{\mathrm{int}}=e\sum_{i\sigma}\mathbf{r}_{i}\cdot\mathbf{E}\left(t\right)c_{i\sigma}^{\dagger}c_{i\sigma},$
(4)

where $\mathbf{E}\left(t\right)=f\left(t\right)E_{0}\hat{\mathbf{e}}\cos\omega
t$ is the electric field strength, with the amplitude $E_{0}$, frequency
$\omega$, polarization $\hat{\mathbf{e}}$ unit vector, and pulse envelope
$f\left(t\right)=\sin^{2}\left(\pi t/\mathcal{T}\right)$. The pulse duration
$\mathcal{T}$ is taken to be $10$ wave cycles: $\mathcal{T}=20\pi/\omega$.
From the Heisenberg equation under the Hartree-Fock approximation one can
obtain evolutionary equations for the single-particle density matrix
$\rho_{ij}^{\left(\sigma\right)}=\left\langle
c_{j\sigma}^{\dagger}c_{i\sigma}\right\rangle$ [65]:

$i\hbar\frac{\partial\rho_{ij}^{\left(\sigma\right)}}{\partial
t}=\sum_{k}\left(\tau_{kj\sigma}\rho_{ik}^{\left(\sigma\right)}-\tau_{ik\sigma}\rho_{kj}^{\left(\sigma\right)}\right)+\left(V_{i\sigma}-V_{j\sigma}\right)\rho_{ij}^{\left(\sigma\right)}$

$+e\mathbf{E}\left(t\right)\left(\mathbf{r}_{i}-\mathbf{r}_{j}\right)\rho_{ij}^{\left(\sigma\right)}-i\hbar\gamma\left(\rho_{ij}^{\left(\sigma\right)}-\rho_{0ij}^{\left(\sigma\right)}\right),$
(5)

where $V_{i\sigma}$ and $\tau_{ij\sigma}$ are defined via density matrix
$\rho_{ij}^{\left(\sigma\right)}$ and its initial value:

$V_{i\sigma}=\sum_{j\alpha}V_{ij}\left(\rho_{jj}^{\left(\alpha\right)}-\rho_{0jj}^{\left(\alpha\right)}\right)+U\left(\rho_{ii}^{\left(\overline{\sigma}\right)}-\rho_{0ii}^{\left(\overline{\sigma}\right)}\right),$
(6)

$\tau_{ij\sigma}=t_{ij}+V_{ij}\left(\rho_{ji}^{\left(\sigma\right)}-\rho_{0ji}^{\left(\sigma\right)}\right).$
(7)

In addition, we assumed that the system relaxes at a rate $\gamma$ to the
equilibrium $\rho_{0ij}^{\left(\sigma\right)}$ distribution. As we see, due to
the mean field modification hopping integrals (7) become non-zero between the
remote nodes, irrespective of the distance.

Figure 1: The top and bottom pannels represent C240 fullerene and C240 with a
monovacancy, respectively. For the brevity for refer to the latter as C239.
Within each row, the following visualizations are presented from left to
right: electron probability density corresponding to the highest energy level
in the valence band on the 3D color mapped molecular structures and the linear
absorption spectra, for Coulomb interaction turned on and off.

##  III Results

Now we discuss full numerical solution of the evolutionary equations for the
single-particle density matrix (5) and to get more physical insight, we study
the question: which effects can be already observed in a linear regime of
interaction. The time propagation of Eq. (5) is performed by the 8-order
Runge-Kutta algorithm. As an initial density matrix we take a fully occupied
valence band and a completely empty conduction band. To study the HHG process
in giant fullerene molecule we evaluate the high-harmonic spectrum by Fourier
transformation of the dipole acceleration,
$\mathbf{a}\left(t\right)=d^{2}\mathbf{d(t)}/dt^{2}$, where the dipole
momentum is defined as
$\mathbf{d}\left(t\right)=e\sum_{i\sigma}\mathbf{r}_{i}\rho_{ii}^{\left(\sigma\right)}\left(t\right)$:

$\mathbf{a}\left(\Omega\right)=\int_{0}^{\mathcal{T}}\mathbf{a}\left(t\right)e^{i\Omega
t}W\left(t\right)dt,$

and $W\left(t\right)$ is the window function to suppress small fluctuations
[88] and to decrease the overall background (noise level) of the harmonic
signal. As a window function we take the pulse envelope $f\left(t\right)$. To
obtain the mean picture which does not depend on the orientation of the
molecule with respect to laser polarization, we take the wave polarization
unity vector as
$\hat{\mathbf{e}}=\left(1/\sqrt{3},1/\sqrt{3},1/\sqrt{3}\right)$.

We begin by examining the effect of vacancy on the states near the Fermi
level. In Fig. 1, electron probability density corresponding to the highest
energy level in the valence band on the 3D color mapped molecular structures
are shown. As is seen from this figure, for a vacancy deffected case we have
state strongly localized around the vacancy. Thus, the presence of single
vacancy also breaks the icosahedral symmetry. To examine intrinsic molecular
transitions, we consider the extreme case of an external electric field that
has the shape of a delta-like impulse in time to excite all electronic
eigenmodes of the systems considered. In this case the relaxation rate is
taken to be very small $\hbar\gamma=0.5\ \mathrm{meV}$ to resolve transitions
as much as possible. The right pannels of Fig. 1 show linear absorption
spectra (in arbitrary units), for Coulomb interaction, turned on and off. The
peaks are intrinsic molecular excitation lines and the area of a particular
peak defines the weight of the oscillator strengths. The effects of the EEI
are similar to those of the fullerene molecule $C_{60}$ molecule [80]. The
Coulomb interaction shift peaks to higher energies, and oscillator strengths
at higher energies have relatively larger weight than in the free electron
case. These effects are due to the fact that the long range Coulomb
interactions (3) give rise to large hopping integrals between the remote nodes
(7) in the Hartree-Fock approximation. For the vacancy defected case the
transitions are overall suppresed compared to intriscic case, although the low
energy transitions are strongly modified. From this figure we also see that
the optical gap in fullerene molecule $\mathrm{C}_{240}$ is approximately
$1.7\ \mathrm{eV}$, which is narrower than that in $\mathrm{C}_{60}$ ($2.8\
\mathrm{eV}$). Notably, in both cases the absorption spectra exhibit many
peaks up to the high energies, suggesting the presence of efficient
multiphoton excitation channels and subsequent high-energy single-photon
transitions. These factors play a significant role in shaping the HHG
spectrum, as we will explore in the following.

Figure 2: The HHG spectra in the strong-field regime in logarithmic scale via
the normalized dipole acceleration Fourier transformation
$a\left(\Omega\right)/a_{0}$ (in arbitrary units) for C240 and for C239. The
laser frequency is $\omega=1.2\ \mathrm{eV}/\hbar$. The spectra are shown for
EEI energy $U=6$ $\ \mathrm{eV}$.

Next, we will study more comprehensive the extreme nonlinear response of giant
fullerene molecule C240 and its vacancy-defected conterpart C239. For all
further calculations, except of Fig. 7, the relaxation rate is taken to be
$\hbar\gamma=0.1\ \mathrm{eV}$. For the convenience, we normalize the dipole
acceleration by the factor $a_{0}=e\overline{\omega}^{2}\overline{d},$ where
$\overline{\omega}=1\ \mathrm{eV}/\hbar$ and $\overline{d}=1\
\mathrm{\mathring{A}}$. The power radiated at the given frequency is
proportional to $\left|\mathbf{a}\left(\Omega\right)\right|^{2}$.

In Fig. 2, we show the typical HHG spectra in the strong field regime
($E_{0}=0.5\ \mathrm{V/A}$) for both molecules. For the C240 molecule, the
presence of inversion symmetry restricts the appearance of only odd harmonics
in the HHG spectrum. In contrast, the introduction of a single vacancy in the
C239 molecule disrupts its icosahedral symmetry, resulting in the prominent
emergence of even-order harmonics with enhanced intensity. Besides, we see
strongly nonlinear picture, where the strength of the $9$th harmonic surpasses
that of the $5$th and $7$th harmonics. Additionally, a distinctive plateau
spanning from the $11$th to the $21$st harmonics exhibits comparable
strengths. Notably, for the C239 molecule, the harmonics near the cutoff
display a slight suppression relative to C240 one. This disparity is
attributed to the differing effectiveness of excitation channels, which favors
enhanced harmonics in the case of C240 molecule (see Fig. 1).

Figure 3: The dependence of the HHG spectra on the wave field frequency is
illustrated for C240 (top) and C239 (bottom) using the normalized dipole
acceleration Fourier transformation, $a\left(\Omega\right)/a_{0}$, plotted on
a logarithmic scale. The wave amplitude is taken to be $E_{0}=0.5\
\mathrm{V/\mathring{A}}$. The relaxation rate is set to $\hbar\gamma=0.1\
\mathrm{eV}$. The EEI energy is U = 6 eV.

Let us now consider the influence of the pump wave frequency on the HHG
process within the energy range of $\hbar\omega=1-2$ $\ \mathrm{eV}$. This
analysis is presented in Fig. 3 that illustrates the frequency-dependent HHG
spectra. Notably, we discern that the position of the cutoff harmonic
$N_{\mathrm{cut}}$ demonstrates a relatively gradual response to changes in
the wave field of frequency $\omega$. Additionally, this cutoff exhibits
distinctive peaks within the mid-frequency range. It’s worth noting that in
atomic HHG processes involving free continua, the cutoff harmonic position
$N_{\mathrm{cut}}\sim\omega^{-3}$ [5]. Furthermore, a noteworthy feature
emerges when considering the C239 molecule: even-order harmonics are
suppressed for higher frequency pump waves. This phenomenon can be attributed
to the fact that with higher frequency pump waves, excitation and
recombination channels predominantly involve highly excited states that still
retain the inversion symmetry. Of particular interest is the plateau region
within the spectra. Here, a pattern of alternating variation in relation to
frequency becomes evident, a hallmark of multiphoton resonant transitions
between the valence and conduction bands. This resonant behavior is further
illuminated by Figs. 4 and 5, where we visualize the dependency of emission
strength for the preplateau harmonics on the pump wave frequency. It is
apparent that these harmonics exhibit resonant behavior. Upon a closer
examination of Fig. 1, we discern that the molecular excitations exhibit peaks
coinciding with these resonant frequencies, providing supplementary evidence
for the multiphoton resonant transitions. For instance, in the case of
molecule C240, the highest peak for the $5$th harmonic emerges at around
$1.3\mathrm{eV}$. This frequency aligns with the local peak at $5\omega\sim
6.5\mathrm{eV}/\hbar$ in Fig. 1, accompanied by multiple excitation channels.
Similarly, considering molecule C239, the peak for the $6$th harmonic is
proximate to $1.18$ $\ \mathrm{eV}$, in acordance with the local peak at
$6\omega\sim 7\mathrm{eV}/\hbar$ in Fig. 1. The peaks displayed in Figs. 4 and
5 correspond with similar peaks in the molecular excitation spectra, as
depicted in Fig. 1.

Figure 4: The dependence of emission strength in the case of C240 for the 3rd,
5th, 7th, and 9th harmonics on the pump wave frequency for the setup of Fig.
3. Figure 5: The dependence of emission strength in the case of C239 for the
2nd, 4th, 6th, and 8th harmonics on the pump wave frequency for the setup of
Fig. 3.

The multiphoton resonance-driven characteristics are further supported by the
evident alteration in the population of energy levels within the valence and
conduction bands, as highlighted in Fig. 6. This figure presents the post-
interaction population distribution of energy levels, demonstrating a marked
departure from the equilibrium distribution. This discrepancy underscores the
substantial impact of multiphoton resonant transitions within the HHG process
of giant fullerene C240 under the influence of intense near-infrared laser
fields.

Continuing our exploration, let us examine the influence of the relaxation
rate on the HHG phenomenon across a span of $\hbar\gamma=0.1-0.2$ $\
\mathrm{eV}$. The corresponding dependencies of the HHG spectra on the
relaxation rate are presented in Fig. 7. It is discernible that HHG exhibits
resistance to relaxation processes, with preplateau harmonics, in particular,
displaying notable robustness.

Figure 6: The residual population of levels for the setup of Fig. 3. Figure 7:
The dependencies of the HHG spectra on the relaxation rate illustrated for
C240 (top) and C239 (bottom). The spectra are shown for EEI energy: U = 6 eV.
The pump wave frequency is $\omega=1.5\ \mathrm{eV}/\hbar$. The wave amplitude
is taken to be $E_{0}=0.5\ \mathrm{V/\mathring{A}}$ The color bar shows the
relaxation rate in $\ \mathrm{eV}/\hbar$.

As have been seen from Fig. 1, the position of molecular excitonic lines and
relative intensities depend on EEI. It is also expected HHG yield change due
to EEI. The latter is shown in Fig. 8, where the HHG spectra in the strong-
field regime for different EEI energies are shown for fullerene C240 molecule.
The similar picture we have for C239 molecule. As is seen, HHG yield strongly
depends on the EEI energy. The inclusion of the Coulomb interaction leads to
two noteworthy characteristics in the HHG spectra: (a) the most prominent
feature is a substantial increase in the HHG signal by several orders of
magnitude near the cutoff regime compared to the case of free quasiparticles.
(b) The cutoff frequency is significantly enhanced. The significant
enhancement in the HHG signal can be explained by the strong modification of
hopping integrals (7) and the resulting level dressing due to the mean field
effect. This observation gains further support from the noticeable prominence
of these features in the case of the giant fullerene C240, in stark contrast
to the behavior observed in C60 molecule [65]. Another notable aspect of the
HHG signals in giant fullerene molecules is their dependence on the size of
the molecule. The HHG signals per particle for C240 and C60 are compared in
Fig. 9. As demonstrated, there is a significant increase in the HHG signal for
C240 molecule, a result also observed for C70 molecule according to previous
studies [65]. This enhancement may be attributed to the density of states,
which is indirectly reflected in Fig. 1 via the absorption spectra. The inset
in Fig. 9 shows the linear absorption spectrum for C60 molecule obtained in
the same way, as in Fig. 1. This figure reveals that C240 molecule has
substantially more transition channels than C60 one.

Figure 8: The comparision of HHG signals for C240 at different EEI energies.
The pump wave frequency is $\omega=1.2\ \mathrm{eV}/\hbar$ and wave amplitude
is 0.5 V/A. The relaxation rate is set to $\hbar\gamma=0.1\ \mathrm{eV}$.
Figure 9: The comparision of HHG signals per particle for C240 and C60. The
pump wave frequency is $\omega=1.1\ \mathrm{eV}/\hbar$ and wave amplitude is
0.5 V/A. The relaxation rate is set to $\hbar\gamma=0.1\ \mathrm{eV}$. The
inset shows the linear absorption spectrum for C60 obtained in the same way as
in Fig. 1.

Finally, note that within the scope of described methodology we have explored
the correlation between the cutoff frequency and the intensity of a pump wave
by analysing the HHG spectra for various intensities. The relationship between
the HHG spectra and the amplitude of the wave field for both giant molecules
is visually represented in Fig. 10. This figure prominently illustrates the
nonlinear connection between the pre-plateau harmonics and the amplitude of
the pump wave. The analysis of obtained results reveals that for high
intensities, the positions of the cutoff harmonics can be adequately described
by scaling with the square root of the field strength amplitude. The solid
lines superimposed on the density plot in Fig. 10, represent envelopes
($\sim\sqrt{E_{0}}$) that determine the positions of the cutoff harmonics.
Notably, it is evident that these envelopes provide a reasonably accurate
approximation for the cutoff harmonics for a large field strengths.

Figure 10: The dependencies of the HHG spectra on the wave field amplitude is
illustrated for C240 (top) and C239 (bottom) using the normalized dipole
acceleration Fourier transformation, $a\left(\Omega\right)/a_{0}$ (color bar),
plotted on a logarithmic scale. The spectra are shown for EEI energy: U = 6
eV. The pump wave frequency is $\omega=1.5\ \mathrm{eV}/\hbar$.

##  IV Conclusion

We have done an extensive exploration of the highly nonlinear optical response
of giant fullerene molecules, with a particular emphasis on C240, which
possesses the characteristic icosahedral point group symmetry often
encountered in such molecular systems. To disclose the complete physical
picture of HHG process on giant fullerene molecules with the mentioned
icosahedral symmetry, we have also investigated a vacancy-defected molecule,
C239. Our investigation employed consistent quantum/analytic and numerical
calculation of the HHG spectra using a mean-field methodology that rigorously
accounts for long-range many-body Coulomb interactions too. Through the
solution of the evolutionary equations governing the single-particle density
matrix we have disclosed resonant effects within the HHG spectra and have
demonstrated the fundamental role of Coulomb interaction in shaping the
intensities of the harmonics. A significant enhancement in HHG yield, as
compared with fullerene molecule C60, has been established. Moreover, our
research has elucidated that the presence of a single vacancy, causing the
breakdown of icosahedral symmetry, stimulates the appearance of pronounced
even-order harmonics. In terms of the dependence of the cutoff harmonics on
the intensity of the wave field, we have established that this relationship
can be approximated with greate accuracy by scaling with the square root of
the amplitude of a pump wave strength.

###### Acknowledgements.

The work was supported by the Science Committee of Republic of Armenia,
project No. 21AG-1C014.

## References

  * Avetissian [2015] H. K. Avetissian,  _Relativistic Nonlinear Electrodynamics: The QED Vacuum and Matter in Super-Strong Radiation Fields_ , Vol. 88 (Springer, 2015).
  * Agostini and DiMauro [2004] P. Agostini and L. F. DiMauro, The physics of attosecond light pulses, Reports on progress in physics 67, 813 (2004).
  * Kohler  _et al._ [2012] M. C. Kohler, T. Pfeifer, K. Z. Hatsagortsyan, and C. H. Keitel, Frontiers of atomic high-harmonic generation, in  _Advances in Atomic, Molecular, and Optical Physics_ , Vol. 61 (Elsevier, 2012) pp. 159–208.
  * Corkum [1993] P. B. Corkum, Plasma perspective on strong field multiphoton ionization, Phys. Rev. Lett. 71, 1994 (1993).
  * Lewenstein  _et al._ [1994] M. Lewenstein, P. Balcou, M. Y. Ivanov, A. LHuillier, and P. B. Corkum, Theory of high-harmonic generation by low-frequency laser fields, Phys. Rev. A 49, 2117 (1994).
  * Corkum and Krausz [2007] P. á. Corkum and F. Krausz, Attosecond science, Nature physics 3, 381 (2007).
  * Krausz and Ivanov [2009] F. Krausz and M. Ivanov, Attosecond physics, Rev. Mod. Phys. 81, 163 (2009).
  * Falcao and Wudl [2007] E. H. Falcao and F. Wudl, Carbon allotropes: beyond graphite and diamond, Journal of Chemical Technology & Biotechnology: International Research in Process, Environmental & Clean Technology 82, 524 (2007).
  * Tiwari  _et al._ [2016] S. K. Tiwari, V. Kumar, A. Huczko, R. Oraon, A. D. Adhikari, and G. Nayak, Magical allotropes of carbon: prospects and applications, Critical Reviews in Solid State and Materials Sciences 41, 257 (2016).
  * Smalley [1997] R. E. Smalley, Discovering the fullerenes, Reviews of Modern Physics 69, 723 (1997).
  * Kroto  _et al._ [1985] H. W. Kroto, J. R. Heath, S. C. O’Brien, R. F. Curl, and R. E. Smalley, C60: Buckminsterfullerene, nature 318, 162 (1985).
  * Kroto and McKay [1988] H. Kroto and K. McKay, The formation of quasi-icosahedral spiral shell carbon particles, Nature 331, 328 (1988).
  * Geim [2009] A. K. Geim, Graphene: status and prospects, Science 324, 1530 (2009).
  * Güçlü  _et al._ [2014] A. D. Güçlü, P. Potasz, M. Korkusinski, P. Hawrylak, _et al._ ,  _Graphene quantum dots_ (Springer, 2014).
  * Gonzalez  _et al._ [1993] J. Gonzalez, F. Guinea, and M. A. Vozmediano, The electronic spectrum of fullerenes from the Dirac equation, Nuclear Physics B 406, 771 (1993).
  * González  _et al._ [1992] J. González, F. Guinea, and M. A. H. Vozmediano, Continuum approximation to fullerene molecules, Physical review letters 69, 172 (1992).
  * York  _et al._ [1994] D. York, J. P. Lu, and W. Yang, Density-functional calculations of the structure and stability of C 240, Physical Review B 49, 8526 (1994).
  * Scuseria [1995] G. E. Scuseria, The equilibrium structures of giant fullerenes: faceted or spherical shape? An ab initio Hartree-Fock study of icosahedral C240 and C540, Chemical physics letters 243, 193 (1995).
  * Scuseria [1996] G. E. Scuseria, Ab initio calculations of fullerenes, Science 271, 942 (1996).
  * Itoh  _et al._ [1996] S. Itoh, P. Ordejón, D. A. Drabold, and R. M. Martin, Structure and energetics of giant fullerenes: An order-N molecular-dynamics study, Physical Review B 53, 2132 (1996).
  * Xu and Scuseria [1996] C. H. Xu and G. E. Scuseria, An O (N) tight-binding study of carbon clusters up to C8640: the geometrical shape of the giant icosahedral fullerenes, Chemical physics letters 262, 219 (1996).
  * Haddon  _et al._ [1997] R. Haddon, G. Scuseria, and R. Smalley, C240—The most chemically inert fullerene?, Chemical physics letters 272, 38 (1997).
  * Heggie  _et al._ [1998] M. Heggie, M. Terrones, B. Eggen, G. Jungnickel, R. Jones, C. Latham, P. Briddon, and H. Terrones, Quantitative density-functional study of nested fullerenes, Physical Review B 57, 13339 (1998).
  * Dunlap and Zope [2006] B. I. Dunlap and R. R. Zope, Efficient quantum-chemical geometry optimization and the structure of large icosahedral fullerenes, Chemical physics letters 422, 451 (2006).
  * Zope  _et al._ [2008] R. R. Zope, T. Baruah, M. R. Pederson, and B. Dunlap, Static dielectric response of icosahedral fullerenes from C 60 to C 2160 characterized by an all-electron density functional theory, Physical Review B 77, 115452 (2008).
  * Calaminici  _et al._ [2009] P. Calaminici, G. Geudtner, and A. M. Koster, First-principle calculations of large fullerenes, Journal of Chemical Theory and Computation 5, 29 (2009).
  * Dunk  _et al._ [2012] P. W. Dunk, N. K. Kaiser, C. L. Hendrickson, J. P. Quinn, C. P. Ewels, Y. Nakanishi, Y. Sasaki, H. Shinohara, A. G. Marshall, and H. W. Kroto, Closed network growth of fullerenes, Nature communications 3, 855 (2012).
  * Martin  _et al._ [2017] J. W. Martin, G. J. McIntosh, R. Arul, R. N. Oosterbeek, M. Kraft, and T. Söhnel, Giant fullerene formation through thermal treatment of fullerene soot, Carbon 125, 132 (2017).
  * Wang  _et al._ [2020] S. Wang, Q. Chang, G. Zhang, F. Li, X. Wang, S. Yang, and S. I. Troyanov, Structural studies of giant empty and endohedral fullerenes, Frontiers in Chemistry 8, 607712 (2020).
  * Ghavanloo  _et al._ [2023] E. Ghavanloo, H. Rafii-Tabar, A. Kausar, G. I. Giannopoulos, and S. A. Fazelzadeh, Experimental and computational physics of fullerenes and their nanocomposites: Synthesis, thermo-mechanical characteristics and nanomedicine applications, Physics Reports 996, 1 (2023).
  * Donnelly  _et al._ [1996] T. D. Donnelly, T. Ditmire, K. Neuman, M. Perry, and R. Falcone, High-order harmonic generation in atom clusters, Physical review letters 76, 2472 (1996).
  * Vozzi  _et al._ [2005] C. Vozzi, M. Nisoli, J. Caumes, G. Sansone, S. Stagira, S. De Silvestri, M. Vecchiocattivi, D. Bassi, M. Pascolini, L. Poletto, _et al._ , Cluster effects in high-order harmonics generated by ultrashort light pulses, Applied Physics Letters 86 (2005).
  * Smirnova  _et al._ [2009] O. Smirnova, Y. Mairesse, S. Patchkovskii, N. Dudovich, D. Villeneuve, P. Corkum, and M. Y. Ivanov, High harmonic interferometry of multi-electron dynamics in molecules, Nature 460, 972 (2009).
  * Mikhailov and Ziegler [2008] S. A. Mikhailov and K. Ziegler, Nonlinear electromagnetic response of graphene: frequency multiplication and the self-consistent-field effects, Journal of Physics: Condensed Matter 20, 384204 (2008).
  * Avetissian  _et al._ [2012] H. Avetissian, A. Avetissian, G. Mkrtchian, and K. V. Sedrakian, Creation of particle-hole superposition states in graphene at multiphoton resonant excitation by laser radiation, Physical Review B 85, 115443 (2012).
  * Avetissian  _et al._ [2013] H. Avetissian, G. Mkrtchian, K. Batrakov, S. Maksimenko, and A. Hoffmann, Multiphoton resonant excitations and high-harmonic generation in bilayer graphene, Physical Review B 88, 165411 (2013).
  * Bowlan  _et al._ [2014] P. Bowlan, E. Martinez-Moreno, K. Reimann, T. Elsaesser, and M. Woerner, Ultrafast terahertz response of multilayer graphene in the nonperturbative regime, Physical Review B 89, 041408(R) (2014).
  * Al-Naib  _et al._ [2014] I. Al-Naib, J. Sipe, and M. M. Dignam, High harmonic generation in undoped graphene: Interplay of inter-and intraband dynamics, Physical Review B 90, 245423 (2014).
  * Chizhova  _et al._ [2016] L. A. Chizhova, F. Libisch, and J. Burgdörfer, Nonlinear response of graphene to a few-cycle terahertz laser pulse: role of doping and disorder, Physical Review B 94, 075412 (2016).
  * Avetissian and Mkrtchian [2016] H. Avetissian and G. Mkrtchian, Coherent nonlinear optical response of graphene in the quantum Hall regime, Physical Review B 94, 045419 (2016).
  * Dimitrovski  _et al._ [2017] D. Dimitrovski, L. B. Madsen, and T. G. Pedersen, High-order harmonic generation from gapped graphene: Perturbative response and transition to nonperturbative regime, Physical Review B 95, 035405 (2017).
  * Avetissian and Mkrtchian [2018] H. Avetissian and G. Mkrtchian, Impact of electron-electron Coulomb interaction on the high harmonic generation process in graphene, Physical Review B 97, 115454 (2018).
  * Sato  _et al._ [2021] S. A. Sato, H. Hirori, Y. Sanari, Y. Kanemitsu, and A. Rubio, High-order harmonic generation in graphene: Nonlinear coupling of intraband and interband transitions, Physical Review B 103, L041408 (2021).
  * Zurrón-Cifuentes  _et al._ [2019] Ó. Zurrón-Cifuentes, R. Boyero-García, C. Hernández-García, A. Picón, and L. Plaja, Optical anisotropy of non-perturbative high-order harmonic generation in gapless graphene, Optics express 27, 7776 (2019).
  * Mrudul and Dixit [2021] M. Mrudul and G. Dixit, High-harmonic generation from monolayer and bilayer graphene, Physical Review B 103, 094308 (2021).
  * Zhang  _et al._ [2021] Y. Zhang, L. Li, J. Li, T. Huang, P. Lan, and P. Lu, Orientation dependence of high-order harmonic generation in graphene, Physical Review A 104, 033110 (2021).
  * Dong  _et al._ [2021] F. Dong, Q. Xia, and J. Liu, Ellipticity of the harmonic emission from graphene irradiated by a linearly polarized laser, Physical Review A 104, 033119 (2021).
  * Avetissian  _et al._ [2022] H. K. Avetissian, G. F. Mkrtchian, and A. Knorr, Efficient high-harmonic generation in graphene with two-color laser field at orthogonal polarization, Phys. Rev. B 105, 195405 (2022).
  * Murakami and Schüler [2022] Y. Murakami and M. Schüler, Doping and gap size dependence of high-harmonic generation in graphene: Importance of consistent formulation of light-matter coupling, Physical Review B 106, 035204 (2022).
  * Tamaya  _et al._ [2023] T. Tamaya, H. Akiyama, and T. Kato, Shear-strain controlled high-harmonic generation in graphene, Physical Review B 107, L081405 (2023).
  * Avchyan  _et al._ [2022a] B. Avchyan, A. Ghazaryan, K. Sargsyan, and K. V. Sedrakian, High harmonic generation in triangular graphene quantum dots, Journal of Experimental and Theoretical Physics 134, 125 (2022a).
  * Avchyan  _et al._ [2022b] B. R. Avchyan, A. G. Ghazaryan, S. S. Israelyan, and K. V. Sedrakian, High harmonic generation with many-body Coulomb interaction in rectangular graphene quantum dots of armchair edge, Journal of Nanophotonics 16, 036001 (2022b).
  * Avchyan  _et al._ [2022c] B. Avchyan, A. Ghazaryan, K. Sargsyan, and K. V. Sedrakian, On Laser-Induced High-Order Wave Mixing and Harmonic Generation in a Graphene Quantum Dot, JETP Letters 116, 428 (2022c).
  * Gnawali  _et al._ [2022] S. Gnawali, R. Ghimire, K. R. Magar, S. J. Hossaini, and V. Apalkov, Ultrafast electron dynamics of graphene quantum dots: High harmonic generation, Phys. Rev. B 106, 075149 (2022).
  * Bauer  _et al._ [2001] D. Bauer, F. Ceccherini, A. Macchi, and F. Cornolti, C 60 in intense femtosecond laser pulses: Nonlinear dipole response and ionization, Physical Review A 64, 063203 (2001).
  * Zhang [2005] G. P. Zhang, Optical high harmonic generation in C 60, Phys. Rev. Lett. 95, 047401 (2005).
  * Zhang and George [2006] G. P. Zhang and T. F. George, Ellipticity dependence of optical harmonic generation in C 60, Phys. Rev. A 74, 023811 (2006).
  * Ciappina  _et al._ [2008] M. Ciappina, A. Becker, and A. Jaroń-Becker, High-order harmonic generation in fullerenes with icosahedral symmetry, Physical Review A 78, 063405 (2008).
  * Ganeev  _et al._ [2009a] R. Ganeev, L. E. Bom, J. Abdul-Hadi, M. Wong, J. Brichta, V. Bhardwaj, and T. Ozaki, Higher-order harmonic generation from fullerene by means of the plasma harmonic method, Physical Review Letters 102, 013903 (2009a).
  * Ganeev  _et al._ [2009b] R. Ganeev, L. E. Bom, M. Wong, J.-P. Brichta, V. Bhardwaj, P. Redkin, and T. Ozaki, High-order harmonic generation from C 60-rich plasma, Physical Review A 80, 043808 (2009b).
  * Redkin and Ganeev [2010] P. Redkin and R. Ganeev, Simulation of resonant high-order harmonic generation in a three-dimensional fullerenelike system by means of a multiconfigurational time-dependent Hartree-Fock approach, Physical Review A 81, 063825 (2010).
  * Topcu  _et al._ [2019] T. Topcu, E. Bleda, and Z. Altun, Drastically enhanced high-order harmonic generation from endofullerenes, Physical Review A 100, 063421 (2019).
  * Ganeev  _et al._ [2013] R. A. Ganeev, C. Hutchison, T. Witting, F. Frank, S. Weber, W. A. Okell, E. Fiordilino, D. Cricchio, F. Persico, A. Zaïr, _et al._ , High-order harmonic generation in fullerenes using few-and multi-cycle pulses of different wavelengths, JOSA B 30, 7 (2013).
  * Zhang and Bai [2020] G. P. Zhang and Y. H. Bai, High-order harmonic generation in solid C 60, Phys. Rev. B 101, 081412(R) (2020).
  * Avetissian  _et al._ [2021] H. K. Avetissian, A. G. Ghazaryan, and G. F. Mkrtchian, High harmonic generation in fullerene molecules, Phys. Rev. B 104, 125436 (2021).
  * Avetissian  _et al._ [2023] H. K. Avetissian, S. Sukiasyan, H. H. Matevosyan, and G. F. Mkrtchian, Disorder-induced effects in high-harmonic generation process in fullerene molecules, Results Phys. 53, 106951 (2023).
  * Orlando  _et al._ [2018] G. Orlando, C.-M. Wang, T.-S. Ho, and S.-I. Chu, High-order harmonic generation in disordered semiconductors, JOSA B 35, 680 (2018).
  * Yu  _et al._ [2019a] C. Yu, K. K. Hansen, and L. B. Madsen, Enhanced high-order harmonic generation in donor-doped band-gap materials, Phys. Rev. A 99, 013435 (2019a).
  * Yu  _et al._ [2019b] C. Yu, K. K. Hansen, and L. B. Madsen, High-order harmonic generation in imperfect crystals, Phys. Rev. A 99, 063408 (2019b).
  * Pattanayak  _et al._ [2020] A. Pattanayak, G. Dixit, _et al._ , Influence of vacancy defects in solid high-order harmonic generation, Physical Review A 101, 013404 (2020).
  * Iravani  _et al._ [2020] H. Iravani, K. K. Hansen, and L. B. Madsen, Effects of vacancies on high-order harmonic generation in a linear chain with band gap, Physical Review Research 2, 013204 (2020).
  * Chinzei and Ikeda [2020] K. Chinzei and T. N. Ikeda, Disorder effects on the origin of high-order harmonic generation in solids, Physical Review Research 2, 013033 (2020).
  * Hansen and Madsen [2022] T. Hansen and L. B. Madsen, Doping effects in high-harmonic generation from correlated systems, Physical Review B 106, 235142 (2022).
  * Xia  _et al._ [2022] C.-L. Xia, J.-Q. Liu, L.-J. Lü, A.-W. Zeng, Z.-L. Li, and X.-B. Bian, Theoretical study of high-order harmonic generation in solutions, Journal of Physics B: Atomic, Molecular and Optical Physics 55, 045401 (2022).
  * Orlando  _et al._ [2022] G. Orlando, M.-I. Lee, and T.-S. Ho, Ellipticity dependence of high-order harmonic generation in disordered semiconductors, Journal of Physics B: Atomic, Molecular and Optical Physics 55, 185601 (2022).
  * Deng  _et al._ [1993] J. P. Deng, D. D. Ju, G. R. Her, C. Y. Mou, C. J. Chen, Y. Y. Lin, and C. C. Han, Odd-numbered fullerene fragment ions from C60 oxides, The Journal of Physical Chemistry 97, 11575 (1993).
  * Banhart [1999] F. Banhart, Irradiation effects in carbon nanostructures, Reports on progress in physics 62, 1181 (1999).
  * Terrones  _et al._ [2000] M. Terrones, H. Terrones, F. Banhart, J.-C. Charlier, and P. Ajayan, Coalescence of single-walled carbon nanotubes, Science 288, 1226 (2000).
  * Martin and Ritchie [1993] R. L. Martin and J. P. Ritchie, Coulomb and exchange interactions in C 60 n-, Phys. Rev. B 48, 4845 (1993).
  * Harigaya and Abe [1994] K. Harigaya and S. Abe, Optical-absorption spectra in fullerenes C 60 and C 70: Effects of Coulomb interactions, lattice fluctuations, and anisotropy, Physical Review B 49, 16746 (1994).
  * [81] M. Yoshida, VRML gallery of Fullerenes (1991), The database is available at http://www. jcrystal. com/steffenweber/gallery/Fullerenes/Fullerenes. html 20.
  * Gilbert [2012] A. Gilbert, IQmol molecular viewer (2012).
  * Ding [2005] F. Ding, Theoretical study of the stability of defects in single-walled carbon nanotubes as a function of their distance from the nanotube end, Phys. Rev. B 72, 245409 (2005).
  * Pereira  _et al._ [2008] V. M. Pereira, J. L. Dos Santos, and A. C. Neto, Modeling disorder in graphene, Phys. Rev. B 77, 115109 (2008).
  * Lee  _et al._ [2005] G.-D. Lee, C. Wang, E. Yoon, N.-M. Hwang, D.-Y. Kim, and K. Ho, Diffusion, coalescence, and reconstruction of vacancy defects in graphene layers, Phys. Rev. Lett. 95, 205501 (2005).
  * Ohno [1964] K. Ohno, Some remarks on the Pariser-Parr-Pople method, Theoretica chimica acta 2, 219 (1964).
  * Harigaya [1998] K. Harigaya, Effects of Coulomb interaction on the nonlinear optical response in, and higher fullerenes, Journal of Physics: Condensed Matter 10, 6845 (1998).
  * Zhang  _et al._ [2018] G. P. Zhang, M. S. Si, M. Murakami, Y. H. Bai, and T. F. George, Generating high-order optical and spin harmonics from ferromagnetic monolayers, Nature Communications 9, 3031 (2018).