11institutetext: Engineering and Public Policy, Carnegie Mellon University, PA
15213, USA 22institutetext: Institute for Software Research, Carnegie Mellon
University, PA 15213, USA 22email<EMAIL_ADDRESS><EMAIL_ADDRESS>
# Climate Change Conspiracy Theories on Social Media

Aman Tyagi 11 0000-0002-6654-0670 Kathleen M. Carley 1122 0000-0002-6356-0238

###### Abstract

One of the critical emerging challenges in climate change communication is the
prevalence of conspiracy theories. This paper discusses some of the major
conspiracy theories related to climate change found in a large Twitter corpus.
We use a state-of-the-art stance detection method to find whether conspiracy
theories are more popular among Disbelievers or Believers of climate change.
We then analyze which conspiracy theory is more popular than the others and
how popularity changes with climate change belief. We find that Disbelievers
of climate change are overwhelmingly responsible for sharing messages with
conspiracy theory-related keywords, and not all conspiracy theories are
equally shared. Lastly, we discuss the implications of our findings for
climate change communication.

###### Keywords:

Climate Change Conspiracy Theories Twitter Conversations Hashtags Label
Propagation Climate Change Communication

##  1 Introduction

There is a virtually 100% consensus among scientists that greenhouse gas
emissions from human activity cause climate change [12]. Despite the
overwhelming evidence, much public discourse shows open skepticism with many
popular contrarian voices [13, 6, 17]. In fact, it is believed that between
20% to 40% of the U.S. population considers climate change as a hoax or do not
believe in its anthropogenic cause [22].

Contrarian voices on climate change can be divided among different categories.
For instance, there is a category of people who argue that climate change is
real but is not caused by human activity. Another example would be people who
believe that climate change is real, but they dispute the anthropogenic cause.
However, the most alarming category are climate change deniers who outright
reject climate science findings or the data as a hoax. Different ideologies
drive most people who describe climate science findings and data as hoax [22].
One such facet of ideology is conspiratorial thinking. Previous studies have
suggested that conspiratorial thinking is associated with beliefs about
climate change. In other words, individuals who believe in conspiracies are
more likely to refute the anthropogenic cause of climate change [22].

Conspiracy theories are “unsubstantiated explanations of events or
circumstances that accuse powerful malevolent groups of plotting in secret for
their own benefit against the common good” [22]. People who believe in
conspiracy theory might want to derive an explanation for any complex
scientific fact from these theories. Conspiracy theories can be interlinked
with each other, although they might not have any logical basis. In this
paper, we would discuss some of the major conspiracy theories in climate
change that are popular on a social media platform. These conspiracy theories
present a significant challenge in removing the false narratives around
climate change. Thus, it becomes essential to analyze these conspiracy
theories.

Previous work on climate change and conspiracy theories suggest that people
believing in conspiracy theories are likely to believe that climate change is
a hoax [22]. That work relied on the manual survey-based collection methods.
Surveys are limited in finding nuanced beliefs and in studying extensive
social network structures. This paper uses extensive Twitter data to link
beliefs about climate change and sharing of conspiracy related text. Thus, the
main research question we answer is, _In climate change discussion, do climate
change Disbelievers share more conspiracy related terms compared to
Believers?_ To answer this research question, we scrape Twitter data for all
the Tweets containing climate change and conspiracy related keywords. We then
use a state-of-the-art stance detection method to find climate change
Believers and Disbelievers 111We define Believers as people who cognitively
accept anthropogenic causes of climate change Disbelievers as those who reject
the same..

Moreover, conspiracies about climate change could be promulgated by bot-like
accounts - automated user accounts - in addition to human actors. These bot-
like accounts can further create confusion on well established climate change
realities. Previous studies have suggested that “bots seek to create false
amplification of contentious issues with the intention to create discord”
[20]. This paper examines whether or not bot-like accounts are more active in
sharing conspiratorial messages in different belief groups.

This paper begins by defining some of the common climate change-related
conspiracy theories §2. Next, we discuss the method used to identify
individual beliefs, keywords used to identify the conspiracy theories, and
method used to find bot-like accounts §3. We present our results in §4. Our
results (§4) suggest that Disbelievers share most conspiracy theory related
Tweets. Conspiracy theory related to chemtrails and geo-engineering is most
popular in our dataset. However, conspiracy theory related to flat earth is
most popular among Believers but rather used as a sarcasm. We also find that
most Disbelievers share only one or two different conspiracy theories with
climate change discussion. Finally, we discuss our findings and their
implications in §5.

##  2 Major Conspiracies about Climate Change

Conspiracy theories evolve with time and do not follow logical arguments. This
section covers the most well-known conspiracy theories related to climate
change and gives brief backgrounds about each. The list was created based on
the author’s findings, and readings of previous work on the same topic [23,
20, 22, 19]222List of example tweets can be found at
https://github.com/amantyag/Climate_Change_Conpiracies/.

  1. 1.

Deep state: Followers of this conspiracy theory agree that there is a hidden
government within the legitimately elected government that controls the state.
Climate change is a hidden agenda of the deep state to further the deep’s
states motives.

  2. 2.

Chem Trails: The condensation trails from the jet engines of an aircraft are
erroneously recognized as consisting of chemical or biological agents. The
theory posits that these trails are responsible for climate change.

  3. 3.

Sunspots: Sunspots are a temporary phenomenon of reduced temperature on the
Sun’s surface [18]. This theory asserts that sunspots and not human activity
are causing climate change.

  4. 4.

Directed Energy Weapon (DEW): A human-made weapon that damages its target by a
highly focussed beam of energy. As per the proponents of this theory, the
usage of DEWs is causing climate change.

  5. 5.

Flat Earth: Advocates of this conspiracy theory do not believe that the earth
is a sphere but rather believe that the earth is a flat disc. Climate is hence
not governed by the standard scientific laws, and climate change is a hoax.

  6. 6.

Geo Engineering: Enthusiasts of this conspiracy theory believe that
governmental experiments cause climate change.

  7. 7.

Unknown Planet: A ninth planet with a vast orbit and unknown to humanity is
causing climate change. The effect of the planet will keep on increasing as it
goes through its perigee.

##  3 Data Collection and Method

In this section, we first describe our data collection in §3.1. Second, in
§3.2, we describe our method to find climate change belief stance and the
keywords used to find conspiracy theory related Tweets.

###  3.1 Data Collection

We collected tweets using Twitter’s standard
API333https://developer.Twitter.com/en/docs/tweets/search/overview/standard
with keywords “Climate Change”, “#ActOnClimate”, “#ClimateChange”. Our dataset
was collected between August 26th, 2017 to September 14th, 2019. Due to server
errors, the collection was paused from April 7th, 2018 to May 21st, 2018, and
again from May 12th, 2019 to May 16th, 2019. We ignore these periods in our
analysis. After deduplicating tweets, our dataset consisted of 38M unique
tweets and retweets from 7M unique users.

###  3.2 Method

This section will first discuss the stance detection method used to identify
climate change Believers and Disbelievers. Then, the keywords used to identify
conspiracy related Tweets.

#### Stance Detection:

Labeling each user as a climate change Believer or a Disbeliever is a non-
trivial task. The broader field of labeling users based on the position the
user takes on a particular topic is called stance mining [16]. We use state-
of-the-art stance mining method which uses weak supervision to find Believers
and Disbelievers [15]. The model uses text signals from Tweets along with
retweet and hashtag network features using a co-training approach with label
propagation [24] and text classification. A set of seed hashtags are provided
as a pro and anti stance signals to the model. The model then labels seed
users based on the usage of these seed hashtags at the end of the tweet
(endtags). The labeled and unlabeled users are then taken as input to the co-
training algorithm. In each step, a combined user-retweet and user-hashtag
network is used to propagate labels to unlabelled users. Concurrently, the
text classifier uses the seed user’s tweets to train an SVM [11] based text
classifier to predict unlabeled users. A common set from text classification
and label propagation of highly confident labels are then used as seed labels
for the next iteration. The final classification is based on the prediction of
the joint model using the combined confidence scores.444 We use the parameter
values as defined in [15] as {$k=5000$, $p=5000$, $\theta^{I}=0.1$,
$\theta^{U}=0.0$, $\theta^{T}=0.7$}. The model has been shown to be above 80%
accurate with multiple datasets.

We select hashtag #ClimateHoax and #ClimateChangeIsNotReal as Disbeliever seed
hashtags and #ClimateChangeIsReal and #SavetheEarth as Believer seed hashtags.
Hashtags ClimateHoax has been shown to be used mostly by Disbelivers [20, 21].
We found similar results on using other Disbeliever hashtags reported in [20,
21]. We use ClimateChangeIsReal and SavetheEarth as Believer hashtags because
of their semantics. Out of the 7M users, we classified 3.1M as disbelievers
and 3.9M as believers 555We randomly sampled 1000 users from each group to
manually validate the results. We label a user as Disbeliever if we find any
Tweet akin to someone who does not believe in climate change or anthropogenic
cause of climate change. Otherwise, we label the user as Believer. We observe
that the average precision from manual validation of 2000 users is 81.2%..

#### Conspiracy Keywords:

We use the following keywords to identify if a Tweet is a conspiracy related
Tweet.

  1. 1.

Deep state: club of rome, clubofrome, clubrome, pizzagate, lizard people,
lizardpeople, illuminati, deepstate, deep state, qanon

  2. 2.

Chem Trails: chemtrail, chem trail

  3. 3.

Sunspots: sunspot

  4. 4.

Directed Energy Weapon: dew, directed energy weapon, directedenergy

  5. 5.

Flat Earth: flat earth, flatearth

  6. 6.

Geo Engineering: geo engineering, geoengineering, weather modification,
weathermodification

  7. 7.

Unknown Planet: planet x, niburu

#### Bot Detection:

We label an account as bot-like or not using CMU’s Bot-Hunter [2, 3]. Bot-
Hunter’s output is a probability measure of bot-like behavior assigned to each
account.

##  4 Results

Climate change Disbelievers share more conspiracy related Tweets than climate
change Believers. We find that the number of Tweets and Retweets shared by
climate change Believers (4830 and 3576) and Disbelievers (31084 and 14369)
are an order of magnitude different. Disbelievers overwhelmingly share Tweets
related to conspiracy theories. Interestingly, for both the groups, conspiracy
theory related Tweets are Tweeted more than Retweeted. This behavior is
contrary to the findings of most studies on Twitter which conclude that users
prefer Retweeting to Tweeting [5]. More Tweeting activity than Retweet
activity suggests that although conspiracy related Tweets can be found in
climate change discussion, not many users are re-sharing the message.

Once we know that Disbelievers are predominantly sharing the conspiracy theory
related Tweets, next, we find which conspiracy theory is most popular. We
break down the Tweets/Retweets with the respective type of conspiracy theory
in figure 1. As expected, Disbelievers are sharing conspiracy theories more
than Believers. The most popular conspiracy theory among Disbelievers is Geo-
engineering and Chemtrails related conspiracy theory. On the other hand,
Believers are sharing Flat Earth conspiracy theory more than other
conspiracies. A manual analysis of 100 randomly selected Tweets show that the
Flat Earth conspiracy theory is used as a sarcastic comment or to make fun of
the other group. We provide further evidence by finding the average sentiment
towards conspiracy related keywords 666To find sentiment towards keywords, we
utilize Netmapper [10] which uses a word-level sentiment computation based on
the average of known valences of surrounding words within a sliding window.
The output values are between -1 and 1, where negative value represents a
negative sentiment, and a positive value represents positive sentiment..
Figure 1 reports the average sentiment in Tweets towards conspiracy theory
related words. Flat Earth conspiracy theory stands out with negative
sentiment, more so when shared by Disbelievers. In other words, irrespective
of beliefs about climate change, the Flat Earth conspiracy theory is viewed
negatively. Interestingly, Believers have a higher positive sentiment towards
ChemTrails and Geo-Engineering conspiracy theories compared to Disbelievers.
We suspect that this could be attributed to Believers explaining the
actualities of these theories. More robust sentiment analysis with a labeled
dataset is needed to draw detailed sentiment level conclusions; such analysis
is out of scope for this work.

Figure 1: Number of unique Tweets and Retweets shared by Disbelievers and
Believers containing different conspiracy theory related keywords defined in
§3.2 (left).Average sentiment score towards the keywords related to different
conspiracy theories. A negative value means a negative sentiment and a
positive value means a positive sentiment towards the conspiracy keywords
(right).

After analyzing the origin of different conspiracy theories, next, we look at
the correlation of different conspiracy theories shared by each user. In table
1, we report the correlation between two different conspiracy theories by
finding the number of times different conspiracy keyword is used by each user.
We find that most conspiracies are highly correlated with each other,
indicating that users who Tweet about one conspiracy also tweet about other
conspiracies. The Chemtrails and Unknown Planet conspiracies are least likely
to be Tweeted by a user who tweets other conspiracy theories. To further gain
insight into the sharing pattern, in figure 2 we report the number of users
sharing unique conspiracy theories. Even on a log scale, we see a steep
decline in the number of Believers and Disbelievers sharing different types of
conspiracy theories.

Next, we find whether or not the same Tweet has more than one conspiracy
discussed. In table 2, we report the correlation between two different
conspiracy theories by finding the number of times different conspiracy
keyword occurs in the same Tweet. We notice in table 2 that there is a weak
negative or close to zero correlation between all the keywords belonging to
different conspiracy theories. Twitter users prefer using conspiracy theory
keywords independent of using other conspiracy theory keywords in a Tweet.
Moreover, conspiracy theories related to Flat Earth and Geo Engineering are
most negatively correlated (-0.338). In other words, Twitter users using
keywords related to Flat Earth do not use keywords related to Geo Engineering
in the same Tweet. Thus, we conclude that most users share one or two types of
conspiracy theory and most Tweets have keywords related to one type of
conspiracy. Moreover, this behavior does not differ from a change in climate
change belief.

Table 1: Correlation matrix of conspiracy theories related keywords used by different users. We find the correlation between two different conspiracy theories by calculating the number of respective keywords used by each user. |  | Deep
---
State
| Chem
---
Trails
Sunspots |  | Directed
---
Energy
Weapons
| Flat
---
Earth
| Geo
---
Engineering
| Planet
---
X
Deep State | 1.000 | 0.360 | 0.713 | 0.953 | 0.958 | 0.892 | 0.161
Chem Trails | 0.360 | 1.000 | 0.195 | 0.333 | 0.321 | 0.361 | 0.041
Sunspots | 0.713 | 0.195 | 1.000 | 0.751 | 0.740 | 0.676 | 0.156
| Directed
---
Energy
Weapons
0.953 | 0.333 | 0.751 | 1.000 | 0.982 | 0.903 | 0.202
Flat Earth | 0.958 | 0.321 | 0.740 | 0.982 | 1.000 | 0.915 | 0.231
| Geo
---
Engineering
0.892 | 0.361 | 0.676 | 0.903 | 0.915 | 1.000 | 0.184
Planet X | 0.161 | 0.041 | 0.156 | 0.202 | 0.231 | 0.184 | 1.000
Table 2: Correlation matrix of conspiracy theories related keywords occurring in a single Tweet. We find the correlation between two different conspiracy theories by calculating the number of respective keywords used in each Tweet. |  | Deep
---
State
| Chem
---
Trails
Sunspots |  | Directed
---
Energy
Weapons
| Flat
---
Earth
| Geo
---
Engineering
| Planet
---
X
Deep State | 1.000 | -0.179 | -0.056 | -0.106 | -0.151 | -0.190 | -0.010
Chem Trails | -0.179 | 1.000 | -0.114 | -0.218 | -0.309 | -0.284 | -0.021
Sunspots | -0.056 | -0.114 | 1.000 | -0.065 | -0.096 | -0.119 | -0.006
| Directed
---
Energy
Weapons
-0.106 | -0.218 | -0.065 | 1.000 | -0.183 | -0.227 | -0.012
Flat Earth | -0.151 | -0.309 | -0.096 | -0.183 | 1.000 | -0.338 | -0.017
| Geo
---
Engineering
-0.190 | -0.284 | -0.119 | -0.227 | -0.338 | 1.000 | -0.022
Planet X | -0.010 | -0.021 | -0.006 | -0.012 | -0.017 | -0.022 | 1.000

Lastly, we look at whether bot-like accounts drive the conspiracy theory
related discussion. We find that even at 0.7 probability cutoff, about a
quarter of all users exhibit bot-like characteristics. We also find that there
is not much difference in the activity between Disbelievers and Believers.
Moreover, we find that most bots ($\sim$88%) share only one type of conspiracy
theory. This conclusion is similar to the results described in figure 2, where
we report the distribution without separating bot-like accounts. Moreover,
bot-like accounts also show a similar pattern with regards to sharing the type
of conspiracy theories. Bot-like accounts showing behavior akin to
Disbelievers share more conspiracy theories related to Geo-engineering and
Chem Trails. On the other hand, bot-like accounts showing behavior akin to
Believers share more conspiracy theories with Flat Earth related keywords.

Figure 2: Distribution of Disbeliever and Believer users (in log scale)
sharing unique conspiracy theories.

##  5 Discussion

Understanding people’s underlying beliefs helps understand the constructs by
which people could be attracted or repelled by different messaging. People
believing in conspiracy theories are more likely to believe that a conspiracy
theory is a possible explanation of climate change [22]. Hence, conspiracy
theories could be used as a potential recruitment tool by Disbeliever
lobbyists. Celebrities and politicians have been vocal about their criticism
of science, even using conspiracy theories as possible explanations for
climate change [22]. These reasons make the study of conspiracy theories in
the climate change context even more relevant.

Conspiracy theories are a means for people to justify the actions of a
powerful entity or a person, mostly when those actions are not relatable [8,
9]. [22] argue that the influence of elites interacting with the masses
predispositions explains conspiracy thinking and why there is a partisan
divide in such thinking. Moreover, President Donald Trump’s election has
further enhanced this effect and could potentially lead to mass radicalization
[1]. Conspiracy belief is thus linked to people’s justification of predisposed
climate change belief. Future research on conspiracy theories warrants these
explanations to be looked at from the lens of psychology and social science.
In this paper, we find that climate change Disbelievers are more likely to
share conspiracy theories. The conspiracy theories range from deep state
conspiracy theory, which portrays climate change as an agenda of individual
actors or deep state to possible explanatory theories such as sunspots and
chemtrails. Future research should look at these theories from the lens of
explanatory or motivated by partisanship.

Conspiracy theories are a real threat to effective climate change messaging.
Climate change messaging should not indiscriminate all conspiracy theories but
tackle the popular ones and alienate the unpopular ones. Climate change
communication research should look to evolve messaging in ways that take into
account different beliefs. Conspiratorial thinking and reasoning to justify
climate change will dampen the global effort to decrease climate change
effects. In this paper, we show that most people sharing conspiracies in the
context of climate change only share one or two types of conspiracies. The
most popular conspiracy theories are related to Chem Trails or Geo-
Engineering. Policymakers should focus on delivering targeted messages to
Disbelievers about the scientific practicalities of these conspiracies.
Moreover, climate change Believers using Flat Earth conspiracy theory to
target Disbelievers or their belief does not help clear scientific facts. Our
results suggest that Flat Earth conspiracy theory is not the most popular
conspiracy theory among Disbelievers.

Previous studies have concluded that Bot-like accounts stir conversations in
differently politically aligned belief groups rather than concentrating on
conversations in one belief group [20, 4]. In this study, we further provide
evidence that Bot-like accounts were similarly active in sharing conspiracy
related messages irrespective of whether they showed activity akin to a
Disbeliever or a Believer. These bot-like accounts aim at widening the divide
between belief groups and pose a danger of creating confusion on scientific
facts [4, 7, 14]. As more and more people consume information via social
media, it becomes imperative for these platforms to identify and remove bot-
like accounts.

To the best of our knowledge, this study is the first attempt to find
conspiracy theories in climate change in a large social media dataset. We find
that some conspiracy theories are more popular and used widely to justify
climate change compared to others. Future psychology and social science
scholarship should divide conspiratorial thinking into different types of
conspiracies. This will help find the underlying constructs and motivations,
knowing which helps target climate change communication messaging.

Besides the demographic representativeness of the data, there are other
limitations in this analysis. First, although we have many tweets about
climate change conspiracies, it does not encompass those interactions that do
not include our collection keywords. Second, we use a proxy of keywords to
classify Tweets as conspiracy related or not. We do not make an effort to find
if sarcasm or negation is used to call out conspiracies; we leave this to
future scholarship. Last, we focused on the conspiracy theories recorded in
media or found during our search. Many more conspiracy theories could be
widespread in the climate change debate. Nevertheless, we believe that we were
able to analyze the main conspiracy theories more widely popular among general
Twitter users.

## References

  * [1] ALLAM, H.: Right-wing embrace of conspiracy is ’mass radicalization,’ experts warn. NPR (2020)
  * [2] Beskow, D., Carley, K.M., Bisgin, H., Hyder, A., Dancy, C., Thomson, R.: Introducing bothunter: A tiered approach to detection and characterizing automated activity on twitter. In: International Conference on Social Computing, Behavioral-Cultural Modeling and Prediction and Behavior Representation in Modeling and Simulation. Springer (2018)
  * [3] Beskow, D.M., Carley, K.M.: You are known by your friends: Leveraging network metrics for bot detection in twitter. In: Open Source Intelligence and Cyber Crime, pp. 53–88. Springer (2020)
  * [4] Bessi, A., Ferrara, E.: Social bots distort the 2016 us presidential election online discussion. First Monday 21(11-7) (2016)
  * [5] Boyd, D., Golder, S., Lotan, G.: Tweet, tweet, retweet: Conversational aspects of retweeting on twitter. In: 2010 43rd Hawaii international conference on system sciences. pp. 1–10. IEEE (2010)
  * [6] Boykoff, M.T., Boykoff, J.M.: Balance as bias: Global warming and the us prestige press. Global environmental change 14(2), 125–136 (2004)
  * [7] Broniatowski, D.A., Jamison, A.M., Qi, S., AlKulaib, L., Chen, T., Benton, A., Quinn, S.C., Dredze, M.: Weaponized health communication: Twitter bots and russian trolls amplify the vaccine debate. American journal of public health 108(10), 1378–1384 (2018)
  * [8] Brotherton, R.: Suspicious minds: Why we believe conspiracy theories. Bloomsbury Publishing (2015)
  * [9] Brotherton, R., French, C.C., Pickering, A.D.: Measuring belief in conspiracy theories: The generic conspiracist beliefs scale. Frontiers in psychology 4, 279 (2013)
  * [10] Carley, L.R., Reminga, J., Carley, K.M.: ORA & Netmapper. In: International Conference on Social Computing, Behavioral-Cultural Modeling and Prediction and Behavior Representation in Modeling and Simulation. Springer (2018)
  * [11] Cortes, C., Vapnik, V.: Support-vector networks. Machine learning 20(3), 273–297 (1995)
  * [12] Doran, P.T., Zimmerman, M.K.: Examining the scientific consensus on climate change. Eos, Transactions American Geophysical Union 90(3), 22–23 (2009)
  * [13] Dunlap, R.E., McCright, A.M., Yarosh, J.H.: The political divide on climate change: Partisan polarization widens in the US. Environment: Science and Policy for Sustainable Development 58(5), 4–23 (2016)
  * [14] Ferrara, E., Varol, O., Davis, C., Menczer, F., Flammini, A.: The rise of social bots. Communications of the ACM 59(7), 96–104 (2016)
  * [15] Kumar, S.: Social Media Analytics for Stance Mining A Multi-Modal Approach with Weak Supervision. Ph.D. thesis, Carnegie Mellon University (2020)
  * [16] Mohammad, S.M., Sobhani, P., Kiritchenko, S.: Stance and sentiment in tweets. ACM Transactions on Internet Technology (TOIT) 17(3), 1–23 (2017)
  * [17] Painter, J., Gavin, N.T.: Climate skepticism in british newspapers, 2007–2011. Environmental Communication 10(4), 432–452 (2016)
  * [18] Ruzmaikin, A.: Origin of sunspots. Space Science Reviews 95(1-2), 43–53 (2001)
  * [19] Tyagi, A.: Challenges in Climate Change Communication on Social Media. Ph.D. thesis, Carnegie Mellon University (2021), https://search.proquest.com/dissertations-theses/challenges-climate-change-communication-on-social/docview/2504529667/se-2?accountid=9902
  * [20] Tyagi, A., Babcock, M., Carley, K.M., Sicker, D.C.: Polarizing tweets on climate change. In: Halil Bisgin, A.H., Dancy, C., , Thomson, R. (eds.) To appear International Conference SBP-BRiMS. Springer (2020)
  * [21] Tyagi, A., Uyheng, J., Carley, K.M.: Affective polarization in online climate change discourse on twitter. In: International Conference on Advances in Social Network Analysis and Mining. IEEE/ACM, Virtual Conference (2020)
  * [22] Uscinski, J.E., Douglas, K., Lewandowsky, S.: Climate change conspiracy theories. In: Oxford Research Encyclopedia of Climate Science. Oxford University Press (2017)
  * [23] Wikipedia contributors: Wikipedia, the free encyclopedia (2020), https://en.wikipedia.org/wiki/List_of_conspiracy_theories, [Online; accessed 22-October-2020]
  * [24] Zhu, X., Ghahramani, Z.: Learning from labeled and unlabeled data with label propagation. Technical Report CMU-CALD-02-107, Carnegie Mellon University (2002)


# Direct stellarator coil optimization for nested magnetic surfaces with
precise quasi-symmetry

Andrew Giuliani<EMAIL_ADDRESS>Courant Institute of Mathematical
Sciences New York University New York, N.Y. 10012-1185  Florian Wechsung
Courant Institute of Mathematical Sciences New York University New York, N.Y.
10012-1185  Antoine Cerfon  Courant Institute of Mathematical Sciences New
York University New York, N.Y. 10012-1185  Matt Landreman  Institute for
Research in Electronics and Applied Physics University of Maryland College
Park, Maryland 20742, USA  Georg Stadler  Courant Institute of Mathematical
Sciences New York University New York, N.Y. 10012-1185

###### Abstract

We present a robust optimization algorithm for the design of electromagnetic
coils that generate vacuum magnetic fields with nested flux surfaces and
precise quasi-symmetry. The method is based on a bilevel optimization problem,
where the outer coil optimization is constrained by a set of inner least-
squares optimization problems whose solutions describe magnetic surfaces. The
outer optimization objective targets coils that generate a field with nested
magnetic surfaces and good quasi-symmetry. The inner optimization problems
identify magnetic surfaces when they exist, and approximate surfaces in the
presence of magnetic islands or chaos. We show that this formulation can be
used to heal islands and chaos, thus producing coils that result in magnetic
fields with precise quasi-symmetry. We show that the method can be initialized
with coils from the traditional two stage coil design process, as well as
coils from a near axis expansion optimization. We present a numerical example
where island chains are healed and quasi-symmetry is optimized up to surfaces
with aspect ratio 6. Another numerical example illustrates that the aspect
ratio of nested flux surfaces with optimized quasi-symmetry can be decreased
from 6 to approximately 4. The last example shows that our approach is robust
and a cold-start using coils from a near-axis expansion optimization.

##  I Introduction

Nested magnetic surfaces are not guaranteed to exist in stellarators, unlike
in tokamaks, and confinement properties of the stellarator can be negatively
affected by the presence of chaos and island chains in the magnetic field.
Furthermore, the calculation of many physics performance metrics that might be
used to optimize stellarators rely on the assumption of nested magnetic
surfaces. MHD stability, neoclassical confinement calculations, and
optimization for quasi-symmetry on surfacesGiuliani  _et al._ (2022a) are
simpler when the field presents nested magnetic surfaces. For all these
reasons, it is crucial to develop stellarator optimization algorithms that are
robust in the presence of chaotic field lines and island chains, and that will
optimize the stellarator to a state with nested magnetic surfaces and precise
quasi-symmetry, i.e., that can target island and chaos healing. By “precise”
quasi-symmetry, we mean that the optimized coil designs generate magnetic
fields that present an accurate approximation of quasi-symmetry where this
property is targeted.

Quasi-symmetry is a favorable property of a magnetic field that ensures the
collisionless trajectories of charged particles up to a certain energy will be
confined Helander (2014). However, it is not completely known how closely this
property can be approximated or realized in non-axisymmetric magnetic fields.
Recently, stellarators with very accurate (though imperfect) quasi-symmetry on
a volume were discovered Landreman and Paul (2022) using a novel stellarator
optimization software framework Landreman  _et al._ (2021).

Island healing was a target during the design of the NCSX coils, where the aim
was to find a plasma equilibrium and coil set such that specific resonant
components of the field from the former and latter cancelled out Hudson  _et
al._ (2002, 2003). In Ref. Zhu  _et al._ , 2019, a spectral analysis was
conducted to manipulate island width size. An optimization-based approach was
recently proposed, which relies on a combined VMEC-SPEC optimization of a
toroidal boundary surface Landreman, Medasani, and Zhu (2021). An initial
configuration with islands and chaos was computed in a stage I optimization
using VMEC. Then, a combined VMEC-SPEC optimization was completed with the
goal of optimizing away the significant island chains present in the initial
state. The VMEC solution was used to optimize for quasi-symmetry, and the SPEC
solution was used to optimize away the magnetic islands by penalizing the
island width via a Greene’s residue computation. The VMEC field was needed
because the quasi-symmetry penalty term relies on the assumption of nested
magnetic surfaces, which are not guaranteed in the SPEC field. The latter
field was needed because island chains cannot be resolved by the VMEC code. It
is also possible to only target island width using SPEC for finite $\beta$
equilibria, without optimizing for quasi-symmetry, which bore promising
resultsBaillod  _et al._ (2022). One disadvantage of the Greene’s residue
approach is that one must know a priori which resonance one would like to
heal. The Greene’s residue approach can be difficult to use if resonances
enter or leave the magnetic configuration at a given iteration of the
optimization, which happens particularly at the start of the optimization
procedure.

Chaos healing also has been the focus of several studies in the past Cary
(1982); Dommaschk (1982); Cary (1984); Hanson and Cary (1984); Cary and Hanson
(1986); Lee  _et al._ (2022). Instead of eliminating localized island chains
corresponding to specific resonances, the idea here is to increase the plasma
volume that is free from generalized stochasticity. This was first done for
analytic vacuum fields, using perturbative analysis Cary (1982); Dommaschk
(1982); Cary (1984). Later, the Greene residue method was applied in an
optimization framework to find coils generating vacuum fields with reduced
chaos Hanson and Cary (1984); Cary and Hanson (1986). These early works did
not consider the simultaneous goal of having magnetic fields with a high level
of quasi-symmetry. More recently, chaos healing was attempted in Ref. Lee  _et
al._ , 2022, where a quadratic flux minimizing (QFM) penalty was added to a
coil optimization objective to favor nested magnetic surfaces further away
from the magnetic axis and optimize away chaotic regions of the field.
However, doing so does not control the quality of quasi-symmetry on low aspect
ratio surfaces.

In this manuscript, we address some of the shortcomings of the previous
approaches. Specifically, the main contributions of this work are: (1) we
outline a robust approach for computing approximations of flux surfaces even
in the presence of island chains and chaos, (2) based on that surface
computation, we present a novel optimization-based approach to island and
chaos healing, (3) our approach optimizes directly the geometry of
electromagnetic coils, rather than a toroidal boundary surface, and (4) the
algorithm promotes precise quasi-symmetry on nested magnetic surfaces.

In our previous work (Ref. Giuliani  _et al._ , 2022a), we outlined a
numerical method to optimize stellarator coils for quasi-axisymmetry under the
assumption that over the course of the optimization, the rotational transform
remained strongly irrational. Despite this strong requirement, we found
precisely quasisymmetric magnetic fields generated by coils. The procedure in
Ref. Giuliani  _et al._ , 2022a relies on the solution to a partial
differential equation (PDE) that can be difficult to solve numerically. In
this work, we propose a least-squares formulation to solve the PDE in a more
robust fashion, resulting in a numerical optimization procedure that can be
used even when nested flux surfaces do not exist. Our approach is formulated
as a bilevel optimization problem, where the outer coil optimization problem
is constrained by the solution to a PDE, solved in a least squares sense in
the inner optimization problem. This is similar in spirit to the work in Ref.
Dudt  _et al._ , 2022, where quasisymmetric stellarators were found by solving
an outer optimization problem constrained by the least squares solution to the
force balance equations computed using the DESC code. We also focus on curl-
free magnetic fields as they are an important first step in the development of
stellarator optimization algorithms. Furthermore, vacuum fields can serve as
suitable initial states in stellarator optimization approaches in which the
plasma pressure is gradually ramped upBoozer (2019).

##  II Computing surfaces

The goal of this section is to outline a robust numerical method for computing
magnetic surfaces in curl-free magnetic fields $\mathbf{B}\in\mathbb{R}^{3}$.
Even though the external magnetic fields that we use here are always generated
by electromagnetic coils, our method is not restricted to fields represented
in this manner. We use a finite-dimensional representation of a toroidal
surface $\bm{\Sigma}_{s}:[0,1)^{2}\rightarrow(x,y,z)$ that satisfies
$n_{\text{fp}}$-rotational symmetry and stellarator symmetry. $n_{\text{fp}}$
stands for the “number of field periods” and indicates the the number of times
the field repeats itself after a full toroidal rotation. The unknowns, or
surface parameters, in this representation are combined into a vector
$\bm{s}\in\mathbb{R}^{n_{s}}$ where $n_{s}$ is the number of parameters that
describe the surface; for details of this representation see Ref. Giuliani
_et al._ , 2022a. Given a vacuum magnetic field $\mathbf{B}$, we seek to
compute a magnetic surface in Boozer coordinates
$\bm{\Sigma}_{s}(\mathbf{s})$, its rotational transform $\iota$, and the
constant $G$, that solve
$\mathbf{r}(\mathbf{s})=[r_{x}(\mathbf{s}),r_{y}(\mathbf{s}),r_{z}(\mathbf{s})]=0$,
whereGiuliani  _et al._ (2022a)

$\mathbf{r}(\mathbf{s}):=G\frac{\mathbf{B}}{\|\mathbf{B}\|}-\|\mathbf{B}\|\left(\frac{\partial\bm{\Sigma}_{s}}{\partial\varphi}+\iota\frac{\partial\bm{\Sigma}_{s}}{\partial\theta}\right).$
(1)

This equation can be derived by equating two different representations of the
magnetic field, which assume that it is curl- and divergence-free Giuliani
_et al._ (2022a); Boozer (2019). Then, with the dual relations, it can be
shown that magnetic surfaces parametrized in Boozer coordinates satisfy (1).
Solutions to this partial differential equation can only be expected when
$\iota$ is strongly irrational. Based on this assumption, we have presented a
pseudospectral approach to solve (1) in Ref. Giuliani  _et al._ , 2022a, and
used that numerical method in an optimization loop to find coils with accurate
quasi-symmetry. The pseudospectral method aimed to find surfaces for which the
residual was exactly zero at a fixed number of collocation points. We called
these surfaces “BoozerExact” surfaces. However, this numerical method can be
brittle when nested flux surfaces do not exist, e.g. in regions with chaotic
field lines and island chains, which occur at places where the rotational
transform is not strongly irrational. The approach that we describe now is
more robust and can be used to determine surfaces even in regions where the
pseudospectral method (BoozerExact) would have difficulty.

Discretizing this partial differential equation, surfaces are computed by
solving the following constrained least squares minimization problem

$\displaystyle\min_{\mathbf{s}}{}$
$\displaystyle\frac{1}{2}\int_{0}^{1}\int_{0}^{1/n_{\text{fp}}}\|\mathbf{r}(\mathbf{s})\|^{2}~{}d\varphi~{}d\theta$
(2) subject to $\displaystyle V(\bm{\Sigma}_{s})-V_{0}=0,$

where
$\mathbf{r}(\mathbf{s}):[0,1)\times[0,1/n_{\text{fp}})\to\mathbb{R}^{3}$,

$\|\mathbf{r}(\mathbf{s})\|^{2}=r_{x}(\mathbf{s})^{2}+r_{y}(\mathbf{s})^{2}+r_{z}(\mathbf{s})^{2},$
(3)

$V_{0}$ is a given target volume, and

$V(\bm{\Sigma}_{s})=\int_{0}^{1}\int_{0}^{1/n_{\text{fp}}}\frac{1}{3}\bm{\Sigma}_{s}\cdot\mathbf{n}~{}d\varphi~{}d\theta,$
(4)

where
$\mathbf{n}=\partial\bm{\Sigma}_{s}/\partial\varphi\times\partial\bm{\Sigma}_{s}/\partial\theta$.
Note that while $\mathbf{n}$ is in the direction of the surface normal, here
is not in general the unit normal. The aim is to solve (1) in a least-squares
sense by minimizing the quadratic residual. The formula for the volume
enclosed by the surface can be derived by recognizing

$V=\int_{D}~{}dx~{}dy~{}dz=\int_{D}\frac{1}{3}\nabla\cdot\vec{r}~{}dx~{}dy~{}dz,$

where $\vec{r}=x\hat{x}+y\hat{y}+z\hat{z}\in D$ and $D$ is the region enclosed
by the toroidal surface $\bm{\Sigma}_{s}$. Applying the divergence theorem to
the right-hand-side, we obtain (4) after substituting
$\vec{r}=\bm{\Sigma}_{s}$.

We use collocation on a tensor grid to approximate the integrals in (2), which
results in the nonlinear least squares problem

$\min_{\mathbf{s}}~{}\frac{1}{6n_{\varphi}n_{\theta}}\|R(\mathbf{s})\|^{2}+\frac{1}{2}w_{v}(V(\bm{\Sigma}_{s})-V_{0})^{2}$
(5)

where

$R(\mathbf{s})=[r_{x}(\mathbf{s})_{1},r_{y}(\mathbf{s})_{1},r_{z}(\mathbf{s})_{1},\ldots,r_{x}(\mathbf{s})_{N_{\text{c}}},r_{y}(\mathbf{s})_{N_{\text{c}}},r_{z}(\mathbf{s})_{N_{\text{c}}}],$

$R(\mathbf{s})\in\mathbb{R}^{3n_{\varphi}n_{\theta}}$ and the indices
correspond to the collocation points. We use a tensor product grid of
$n_{\varphi}$ and $n_{\theta}$ collocation points in the $\varphi$ and
$\theta$ directions, respectively. We have experimented with two grids of
quadrature points of the form
$(\varphi_{i},\theta_{j})=(i\Delta\varphi,j\Delta\theta)$, for
$i=0,1,\ldots,n_{\varphi}-1$, $j=0,1,\ldots,n_{\theta}-1$. Rule 1 is on a
full-period $(\varphi,\theta)\in[0,1/n_{\text{fp}})\times[0,1)$,
$\Delta\varphi=(1/n_{\text{fp}})/n_{\varphi}$ and $\Delta\theta=1/n_{\theta}$,
$n_{\varphi}=2n_{\text{tor}}+1$, $n_{\theta}=2n_{\text{pol}}+1$. This rule is
spectrally accurate. Rule 2 is on a half-period
$(\varphi,\theta)\in[0,1/2n_{\text{fp}})\times[0,1)$ with
$\Delta\varphi=(1/2n_{\text{fp}})/n_{\varphi}$, $\Delta\theta=1/n_{\theta}$,
$n_{\varphi}=n_{\text{tor}}+1$, $n_{\theta}=2n_{\text{pol}}+1$. This rule
exploits stellarator symmetry, using half the number of points as rule 1 and
is not spectrally accurate. A comparison of the two rules follows in Section
II.1.

A related concept is that of quadratic-flux-minimizing (QFM) surfaces Dewar,
Hudson, and Price (1994), which are surfaces that minimize

$\int_{0}^{1}\int_{0}^{1/n_{\text{fp}}}(\mathbf{B}\cdot\mathbf{n})^{2}~{}d\theta
d\varphi,$

without constraints on the angles that parametrize the surface. QFM and
BoozerLS surfaces are related to one another in regimes where nested flux
surfaces exist. In infinite dimensions, BoozerLS surfaces are also QFM
surfaces, with the added requirement that the BoozerLS surface is
parameterized in Boozer coordinates. In finite dimensions, numerical
experiments show that BoozerLS surfaces do indeed approximate magnetic
surfaces and we expect both BoozerLS and QFM to be close to each other. The
existence of QFM surfaces in more general contexts is delicate: when the
quadratic flux is unweighted, it is shown in Ref. Dewar, Hudson, and Price,
1994 that only true flux surfaces extremize the QFM functional. An analogous
conclusion may be true for BoozerLS surfaces, but we have not attempted to
show this here.

The first-order optimality condition of (5) is

$\mathbf{g}(\mathbf{s}):=J(\mathbf{s})^{T}R(\mathbf{s})+w_{v}(V(\mathbf{s})-V_{0})\frac{\partial
V}{\partial\mathbf{s}}=\mathbf{0}$ (6)

where $\mathbf{g}\in\mathbb{R}^{n_{s}}$, $J=\frac{\partial
R}{\partial\mathbf{s}}\in\mathbb{R}^{(3n_{\varphi}n_{\theta})\times n_{s}}$ is
the Jacobian of $R(\cdot)$ and

$V(\mathbf{s}):=\frac{1}{3n_{\varphi}n_{\theta}}\sum_{i=1}^{3n_{\varphi}n_{\theta}}(\bm{\Sigma}_{s})_{i}\cdot\mathbf{n}_{i}.$
(7)

is the discretized surface volume.

In (Ref. Giuliani  _et al._ , 2022a), we showed numerically that computing
BoozerLS surfaces is robust, even in the presence of islands and chaos.
Framing the surface computation as a least-squares optimization problem allows
us to use line search algorithms that track progress of the solution
algorithm, and to prevent step sizes that are too large. Thanks to the
optimization framework which defines the BoozerLS surfaces, we are free to
introduce regularizations that prevent the numerically computed surfaces from
self-intersecting.

The Newton step of (6) for the increment $\delta\mathbf{s}$ is

$\displaystyle\biggl{[}J(\mathbf{s})^{T}J(\mathbf{s})$
$\displaystyle+w_{v}\frac{\partial V}{\partial\mathbf{s}}^{T}\frac{\partial
V}{\partial\mathbf{s}}$ (8)
$\displaystyle+\sum_{i=1}^{3n_{\varphi}n_{\theta}}R_{i}(\mathbf{s})\frac{\partial^{2}R_{i}}{\partial\mathbf{s}^{2}}(\mathbf{s})$
$\displaystyle+w_{v}(V(\mathbf{s})-V_{0})\frac{\partial^{2}V}{\partial\mathbf{s}^{2}}\biggr{]}\delta\mathbf{s}=-\mathbf{g}(\mathbf{s})$

where the Hessian matrix on the left-hand side multiplying
$\delta\mathbf{s}\in\mathbb{R}^{n_{s}}$ is denoted by
$H\in\mathbb{R}^{n_{s}\times n_{s}}$. We initially use the
Broyden–Fletcher–Goldfarb–Shanno (BFGS) algorithm to solve (5). Then, using
this solution as an initial guess, we use Newton’s method on (8) to reduce the
gradient of the nonlinear residual further. We observe that while the Hessian
is not absolutely necessary for the computation of the least square surfaces,
its availability is crucial for the computation of gradients in the outer coil
optimization problem we consider in Section III.2, when the least-squares
optimality condition is used as constraint within the coil optimization
problem. We will highlight this point again in that section. We also note that
in contrast to Ref. Giuliani  _et al._ , 2022a, the residual in (3) is divided
by $\|\mathbf{B}\|$. This is to prevent $\mathbf{r}$ from scaling with the
coil currents.

The approach described in this section shows how to compute a magnetic surface
that encloses a user-defined volume $V_{0}$. The coil optimization studies in
Section III.2 use this numerical method to compute multiple magnetic surfaces,
each with a different target volume.

###  II.1 Convergence study

In this section, we present a convergence test of the BoozerLS formulations
using two different quadrature rules to approximate the integral in (2),
called rules 1 and 2. We compare this new least squares formulation to the one
in Ref. Giuliani  _et al._ , 2022a i.e., when the number of unknowns coincides
with the number of collocation points, which we refer to as BoozerExact. We
compute the innermost surface used in the coil optimization problem (section
IV) using the BoozerExact and constrained BoozerLS formulations and vary the
number of modes, $m_{\text{pol}},n_{\text{tor}}$, in the surface
representation and the number of collocation points. Note that
$m_{\text{pol}}$ and $n_{\text{tor}}$ are the number of Fourier modes used to
represent the surface in the $\theta$ and $\varphi$ directions, respectively.
In addition to computing the unknowns that describe the geometry of the
surface, we also compute its rotational transform $\iota$ as well as $G$. For
the configuration here, we know the exact value of $G=\mu_{0}\sum_{k}I_{k}$,
where $I_{k}$ is the current in the $k$th coil. Thus, we can plot the
convergence of the numerically computed $G$ to the true value, shown in figure
1. We find that the accuracy of the BoozerLS formulation depends on the choice
of quadrature rule used to approximate the quadratic residual in (6). If a
spectrally accurate quadrature rule is used (rule 1), then BoozerLS is as
accurate as the BoozerExact formulation or better. If a non-spectrally
accurate quadratule is used (rule 2), then BoozerLS is orders of magnitude
less accurate than the BoozerExact formulation.

Since the the BoozerLS formulation is most useful in regions without nested
flux surfaces, where (1) does not even have well-defined solutions, we use
rule 2 in what follows as it requires fewer points than rule 1.

Figure 1: Convergence of $G$ with respect to number of surface modes on the
innermost surface of the initial configuration in section IV.1 and IV.2.
Exponential convergence results in linear error curves on a $\log$-linear plot
as observed here. Due to the spectral accuracy of the quadrature rule, we find
that BoozerLS surfaces computed with rule 1 are more accurate than rule 2.

###  II.2 Surface regularization

Numerical evidence suggests that surfaces determined by (2) converge
exponentially as a function of the number of Fourier modes that describe the
surface in regions with nested flux surfaces. However, we have not examined
what happens in regions with chaos and islands. As we increase the Fourier
resolution, we observe that least squares surfaces may present self-
intersections. One possible surface regularization to prevent the formation of
self-intersections is

$\frac{1}{2}w_{v}(A(\mathbf{s})-A_{0})^{2},$

where $A_{0}$ is 0.9 times the surface area on the unregularized surface, and
$w_{v}>0$ is a weighting parameter for the regularization term. The goal of
this regularization is to prevent the area of the surface from becoming too
large for a given toroidal volume, thereby avoiding self-intersections. The
effect of this regularization term is illustrated in Figure 2, where we
compute a surface that passes through the $\iota=2/5$ island chain, for
$m_{\text{pol}},n_{\text{tor}}=15,16$. We observe that increasing the surface
complexity from $m_{\text{pol}},n_{\text{tor}}=15$ to
$m_{\text{pol}},n_{\text{tor}}=16$, results in surfaces that self-intersect in
the neighborhood of the $X$-points.

Figure 2: Cross sections for $m_{\text{pol}},n_{\text{tor}}=15,16$ surfaces,
plotted with solid and dashed lines, respectively. The red and blue cross
sections correspond to the unregularized and regularized surfaces,
respectively in the initial configuration of section IV.1 and IV.2. The
unregularized surface present a sharp cusp and self-intersections at the
$X$-points, while the regularized surfaces do not. Moreover, the cross
sections of the regularized surfaces overlap closely and cannot be
distinguished visually.

We also compute the spectrum of the Hessian for the
$m_{\text{pol}},n_{\text{tor}}=16$ surface (Figure 3).

Figure 3:  Spectrum of the Hessian at the minimizer of (5) when
$m_{\text{pol}},n_{\text{tor}}=16$. The are many small eigenvalues and as a
result, there are many directions in which the surface can be perturbed
without substantially increasing the magnitude of the residual
$\|R(\mathbf{s})\|^{2}$.

We find that before the regularization is added, there are many small
eigenvalues and as a result, there is much freedom to add design requirements
on the surface without affecting the magnitude of the Boozer residual. This is
confirmed by our numerical tests, as we observe that the regularization
negligibly affects the Boozer residual penalty (Figure 4). The spectrum
changes negligibly after adding the area penalty term, meaning that the door
is open for additional regularizations to be included if needed. The
regularization introduced here is only a heuristic and there are no guarantees
that self-intersections will always be prevented. We find that this solution
is sufficient for our purposes, and useful at the start of a coil optimization
procedure or if a BoozerLS surface is in the neighborhood of a low-order
rational. An alternative regularization is to penalize the surface curvatures,
but we have not pursued this here for simplicity.

Figure 4:  The convergence of the residual $\|R(\mathbf{s})\|^{2}$ with
increasing surface complexity $m_{\text{pol}},n_{\text{tor}}$ with and without
surface area regularization. As predicted by the eigenvalue analysis in Figure
3, the residual on the regularized surface still behaves like the
unregularized one.

In what follows, the surface area regularization is used in section IV.1 and
IV.2 to ensure robustness of coil optimization problems as the computed
BoozerLS surfaces have a high number of Fourier modes
($m_{\text{pol}},n_{\text{tor}}=15$). Section IV.3 does not require the
surface area regularization for robustness, as BoozerLS surfaces with a lower
number of Fourier modes ($m_{\text{pol}},n_{\text{tor}}=10$) are computed.

##  III Coil optimization for quasi-symmetry on surfaces

In this section, we show how to evaluate the quality of quasi-symmetry on
surfaces, and we formulate a coil optimization problem that targets quasi-
symmetry on those surfaces. Finally, we discuss the efficient computation of
gradients for this optimization problem. In the previous section, we were
concerned with developing a numerical method to approximate a single magnetic
surface that encloses a user-provided volume. In what follows, we use the
numerical method to compute $N_{s}$ magnetic surfaces that enclose different
volumes, and use those surfaces in a coil optimization procedure to control
quasi-symmetry in the generated magnetic field.

###  III.1 Measuring quasi-symmetry on surfaces

Once a surface is known, we compute the deviation from quasi-symmetry on that
surface following the approach from (Ref. Giuliani  _et al._ , 2022a). In this
work, we focus only on quasi-axisymmetry (QA), however, there is no
fundamental restriction preventing the extension of this work to other types
of quasi-symmetry. First, the field strength $B(\varphi,\theta)$ on the
surface is decomposed into a quasisymmetric and non-quasisymmetric component

$B(\varphi,\theta)=B_{\text{non-QA}}(\varphi,\theta)+B_{\text{QA}}(\theta).$

The quasisymmetric field strength is computed using a least squares projection

$B_{\text{QA}}(\theta)=\frac{\int^{1/n_{\text{fp}}}_{0}B(\bm{\Sigma}_{s}(\varphi,\theta))~{}\|\frac{\partial\bm{\Sigma}_{s}}{\partial\varphi}\times\frac{\partial\bm{\Sigma}_{s}}{\partial\theta}\|~{}d\varphi}{\int^{1/n_{\text{fp}}}_{0}~{}\|\frac{\partial\bm{\Sigma}_{s}}{\partial\varphi}\times\frac{\partial\bm{\Sigma}_{s}}{\partial\theta}\|~{}d\varphi},$

and is the closest quasisymmetric field strength to $B(\varphi,\theta)$ when
measured in the norm $(\int_{\bm{\Sigma}_{s}}f^{2}~{}dS)^{1/2}$. The objective
of the optimization problem presented in the next section uses this measure of
non-quasi-symmetry to design stellarator coils.

###  III.2 The coil optimization problem

Next, we formulate and solve an optimization problem for magnetic coils such
that the induced magnetic surfaces have good quasi-symmetry. Where magnetic
surfaces do not exist, our least squares framework will still provide a
surface that can be used in the optimization procedure. This is in contrast to
our previous workGiuliani  _et al._ (2022a), which would have had issues in
regions in which the magnetic field does not possess nested flux surfaces. At
the end of the optimization algorithm, islands and generalized chaos will be
healed. We search for $N_{c}$ independent modular coils to which stellarator
and rotational symmetries are applied such that the stellarator is made up of
$2n_{\text{fp}}N_{c}$ modular coils. Each coil is parametrized with a Fourier
representation and a current with $n_{c}$ degrees of freedom as in (Ref.
Giuliani  _et al._ , 2022b). All coil degrees of freedom are summarized into a
vector $\mathbf{c}\in\mathbb{R}^{N_{c}n_{c}}$. We target quasi-symmetry of the
induced magnetic field on $N_{s}$ surfaces, which are characterized as
minimizers of (5) with the first-order necessary conditions

$\mathbf{g}_{k}(\mathbf{s}_{k})=0,k=1,\ldots,N_{s}$

as defined in (6). Each surface satisfies the optimality condition, but has a
different target volume. In this way, we are able to extend the single-surface
method described in Section II to multiple surfaces. The corresponding surface
is denoted by $\bm{\Sigma}_{s,k}$. Here and in the following, the index $k$
corresponds to the $k$th surface. All surface parameters are combined into a
vector $\mathbf{s}:=(\mathbf{s}_{1},\ldots,\mathbf{s}_{N_{s}})$. The objective
that we minimize is the sum of the average (normalized) non quasi-axisymmetry
and the Boozer residuals on the surfaces $\bm{\Sigma}_{s,k}$:

$\displaystyle\hat{f}(\mathbf{c},\mathbf{s})$
$\displaystyle:=\frac{1}{N_{s}}\sum^{N_{s}}_{k=1}\biggl{\\{}\frac{\int_{\bm{\Sigma}_{s,k}}B_{\text{non-
QA}}(\mathbf{c},\mathbf{s}_{k})^{2}~{}dS}{\int_{\bm{\Sigma}_{s,k}}B_{\text{QA}}(\mathbf{c},\mathbf{s}_{k})^{2}~{}dS}$
(9)
$\displaystyle+\frac{1}{2}w_{r}\int_{0}^{1}\int_{0}^{1/n_{\text{fp}}}\|\mathbf{r}_{k}(\mathbf{s})\|^{2}~{}d\varphi~{}d\theta\biggr{\\}},$

where $w_{r}>0$ is a weighting parameter for the Boozer residual. If the
Boozer residual is not in this objective, i.e. $w_{r}=0$, then there is
nothing preventing the accuracy of the least squares surface from degrading as
the coils are optimized. Therefore, we include it to ensure that the least
squares surface residual remains small, which is particularly important in the
presence of islands and chaos, i.e., when no magnetic surfaces exist and the
surface found with (5) is an approximate surface. This term is not necessary
for BoozerExact surfaces, when the same number of surface parameters and
collocation points is used as in (Ref. Giuliani  _et al._ , 2022a).

This residual term can successfully detect regions of the magnetic field
without nested flux surfacesGiuliani  _et al._ (2022a). When the residual term
is small, then the BoozerLS surfaces accurately solve the PDE. When it is
large, then it is more likely that nested surfaces do not exist. The addition
of the residual term in the outer optimization problem will favor coils that
produce nested magnetic surfaces, which solve (1) accurately. This approach is
similar to the one taken in Ref. Lee  _et al._ , 2022, where the quadratic
flux minimizing (QFM) surface penalty was used to recover nested fluxed
surfaces away from the magnetic axis without controlling for the quality of
quasi-symmetry on those surfaces. The technique that we adopt here differs in
that not only can we favor nested flux surfaces with lower aspect ratio, but
we can also directly optimize for quasi-symmetry on those surfaces.

Our goal is to find a set of coils that solves the following optimization
problem:

$\displaystyle\min_{\mathbf{c}\in\mathbb{R}^{N_{c}n_{c}-1},\>\mathbf{s}\in\mathbb{R}^{N_{\\!s}n_{s}}}$
$\displaystyle~{}\hat{f}(\mathbf{c},\mathbf{s})$ (10a)
$\displaystyle\text{subject to }\mathbf{g}_{k}(\mathbf{s}_{k})=0$ for
$k=1,\ldots,N_{s}$, (10b) $\displaystyle c(\iota_{1},\ldots,\iota_{N_{s}})$
$\displaystyle=0,$ (10c) $\displaystyle R_{\text{major}}$
$\displaystyle=R_{0},$ (10d) $\displaystyle\sum_{i=1}^{N_{c}}L^{(i)}_{c}$
$\displaystyle\leq L_{\max},$ (10e) $\displaystyle\kappa_{i}$
$\displaystyle\leq\kappa_{\max},~{}i=1,\ldots,N_{c},$ (10f)
$\displaystyle\frac{1}{L^{(i)}_{c}}\int_{\bm{\Gamma}^{(i)}}\kappa_{i}^{2}~{}dl$
$\displaystyle\leq\kappa_{\mathrm{msc}},~{}i=1,\ldots,N_{c},$ (10g)
$\displaystyle\|\bm{\Gamma}^{(i)}-\bm{\Gamma}^{(j)}\|$ $\displaystyle\geq
d_{\min}~{}\text{ for }i\neq j,$ (10h)
$\displaystyle\|\bm{\Gamma}^{\prime(i)}\|-L^{(i)}$ $\displaystyle=0~{}\text{
for }i=1,\ldots,N_{c},$ (10i)

where $c(\iota_{1},\ldots,\iota_{N_{s}})$ is an equality constraint on the
rotational transform profile that will be specified in the numerical examples.
The optimization is subject to the surface constraints (10b)–(10d) and the
coil constraints (10e)–(10i). Note that in (10b), we use the first-order
optimality conditions (6) rather than the least squares minimization problem
(5) to define surfaces. This is a standard approach to make bilevel
optimization problems, i.e., optimization problems where the constraint is
itself an optimization problem, computationally tractable. Next, (10c)
enforces that the rotational transform on a specific surface or that the
average rotational transform is $\overline{\iota}$; other constraints on the
rotational transform can easily be incorporated. The constraint (10d) fixes
the major radius on the innermost surface to a given $R_{0}$ and prevents the
length scale of the stellarator from changing. (10e) prevents the sum of the
independent modular coil lengths $\sum_{i=1}^{N_{c}}L^{(i)}_{c}(\mathbf{c})$
from exceeding a given value $L_{\max}>0$. (10f) and (10g), respectively,
prevent the curvature and mean squared curvature on each coil from exceeding
the values $\kappa_{\max}$ and $\kappa_{\mathrm{msc}}$. (10h) ensures that the
coils stay at least $d_{\min}>0$ away from one another and (10i) enforces that
the coil parameterization has a uniform arclength. Finally, in order to
prevent the coil currents from approaching zero, the current in the first coil
is fixed to a nonzero value over the course of the coil optimization, thus the
dimensions of the outer and inner optimization problems are respectively
$N_{c}n_{c}-1$ and $N_{s}n_{s}$.

By virtue of the finite dimensional Fourier representation of the surfaces, we
are restricted to magnetic surfaces that have some smoothness. Thus, our
optimized configurations will not likely present divertor surface shapes, even
though it is known that they can be efficiently produced Landreman and Boozer
(2016).

###  III.3 Handling of constraints and gradient computation

The (in)equality constraints (10c)-(10h) are enforced by adding quadratic
penalties to the objective. For example, to constrain the rotational transform
on an inner surface to $2/5$, a quadratic penalty of the form
$(\iota-2/5)^{2}$ is added to the objective, where $\iota$ is the rotational
transform on an interior surface. In practice, the weights that multiply these
quadratic penalties are only increased if the constraint is violated by more
than 0.1%. Next, we discuss the computation of the discretely exact gradient
of the reduced objective
$f(\mathbf{c}):=\hat{f}(\mathbf{c},\mathbf{s}(\mathbf{c}))$, i.e., the
objective in which we consider the surfaces as a function of the coils through
the magnetic field induced by the coils. We use an adjoint method to compute
the gradient $\nabla_{\mathbf{c}}f$ efficiently as follows:

$\displaystyle\nabla_{\mathbf{c}}f$
$\displaystyle=\frac{\partial\hat{f}}{\partial\mathbf{c}}-\sum_{k}^{N_{s}}\lambda_{k}^{T}\frac{\partial\mathbf{g}_{k}}{\partial\mathbf{c}},$
(11a) $\displaystyle H_{k}^{T}\lambda_{k}$
$\displaystyle=\frac{\partial\hat{f}}{\partial\mathbf{s}_{k}}~{}\text{ for
}k=1,\ldots,N_{s}.$ (11b)

Note that solving the $N_{s}$ adjoint systems (11b)—one for each
surface—requires the Hessian matrix from the Newton system (8). This is a
consequence of the use of a least-squares formulation for the surface
computation. Note that we typically finish the surface computation with a few
Newton iterations, and thus the Hessian matrices $H_{k}$ needed in (11) are
already available.

###  III.4 Computational aspects

Many of the tools presented here are implemented in the SIMSOPT software
packageLandreman  _et al._ (2021), which is a suite for stellarator design and
optimization.

The bilevel optimization problem described above can be computationally
expensive. To illustrate this, consider first the inner optimization problems,
which must be solved at each iteration of the outer optimization problem. If
$m_{\text{pol}},n_{\text{tor}}=N$, then there are $\mathcal{O}(N^{2})$ surface
degrees of freedom that are determined by the inner optimization problem.
Since $n_{\varphi},n_{\theta}=\mathcal{O}(N)$, the number of residuals in the
least squares objective scales like $\mathcal{O}(N^{2})$. As a result, the
computational work to evaluate the gradient and Hessian of the Boozer residual
in (2) scales like $\mathcal{O}(N^{4})$ and $\mathcal{O}(N^{6})$,
respectively. These scalings are due to the use of a (globally defined)
Fourier basis in the surface representation. For example, increasing $N=10$ to
$N=15$ will make the gradient and Hessian evaluation approximately 5 and 10
times more computationally expensive, respectively. For these surface
computations, we typically first use BFGS to robustly obtain a good
approximate BoozerLS surface, and then improve the obtained solution using
very few iterations of Newton’s method. As a result, the surface computation
is largely dominated by the BFGS iterations. At the expense of accuracy, some
of the ill-scaling of the computational work could be mitigated by using more
compact basis functions, e.g., finite element polynomial bases, which are
generally not globally defined. We also compared the performance of L-BFGS and
BFGS algorithms for the surface computation, and, while each iteration of
L-BFGS was faster, it typically required more iterations to reach the same
accuracy as BFGS.

Consider now the outer coil optimization problem, which we solve using again
the BFGS algorithm. The computational work to evaluate the value of the
objective and its gradient is dominated by the inner BFGS optimization and we
observe that a single outer iteration takes under a minute when $N=10$ and
over 5 minutes when $N=15$ on Intel Xeon Platinum 8268 processors, which
agrees with the rough scaling computation in the previous paragraph. This is
in contrast to the local optimization algorithm presented in Ref. Giuliani
_et al._ , 2022a, where each iteration of the outer coil optimization
algorithm took on the order of a few seconds when $N=10$.

Much effort has been made to accelerate the optimization algorithm. MPI
parallelism is used when quasi-symmetry is optimized on multiple surfaces,
where each surface computation is completed on separate MPI ranks. On each
rank, the BoozerLS objective, gradient, and Hessian are evaluated using
multiple cores with OpenMP in addition to SIMD parallelism on a given core.

##  IV Experiments

In the following numerical experiments, we aim to demonstrate the robustness
of the BoozerLS formulation and its ability to optimize for quasi-symmetry,
even in the presence of islands and chaos. The general optimization procedure
we follow is to initially use BoozerLS surfaces everywhere. Then, once the
optimizer reaches the neighborhood of a local minimizer, surfaces in the
neighborhood of low-order rationals continue to use the regularized BoozerLS
formulation, while surfaces away from low-order rationals use the BoozerExact
formulation. The magnetic field that we design has two-fold rotational
symmetry ($n_{\text{fp}}=2$), stellarator symmetry, and a major radius
$R_{0}=$1\text{\,}\mathrm{m}$$. For all examples that follow, we set the coil
design requirements to be $d_{\min}=$0.1\text{\,}\mathrm{m}$$,
$\kappa_{\max}=$5\text{\,}\mathrm{m}$^{-1}$,
$\kappa_{\text{msc}}=$5\text{\,}\mathrm{m}$^{-2}$, and
$L_{\max}=$18\text{\,}\mathrm{m}$$. We constrain the rotational transform
profile to pass through the low-order rational $\iota=2/5$ and optimize for QA
up to surfaces with aspect ratio 6 and 4 in sections IV.1 and IV.2,
respectively. In section IV.3, we constrain the average rotational transform
in the volume confined by the outermost surface of aspect ratio 6 to be
$0.42$. Depending on the example, we optimize using eight or nine surfaces.
Using fewer surfaces might have resulted in comparable levels of quasi-
symmetryGiuliani  _et al._ (2022a), but we did not examine this question here.

In all examples, the $x,y,z$ coordinates of the independent modular coils are
described using 16 Fourier modes resulting in $3(2\times 16+1)=99$ geometric
degrees of freedom and an associated current per coil, i.e. the number of
degrees of freedom per coil is $n_{c}=100$. We design stellarators composed of
four independent coils that have $4\times(99+1)=400$ degrees of freedom. The
outer coil optimization problem has $399$ dimensions however, since we fix the
current on the first coil to be constant, thereby preventing the currents of
the stellarator from approaching zero. A complete set of 16 coils is obtained
by applying twofold discrete rotational symmetry and stellarator symmetry. For
example IV.1 and IV.2, we consider surfaces with
$m_{\text{pol}},n_{\text{tor}}=15$, which results in a 1,439 dimensional inner
optimization problem for each surface. For example IV.3, we consider surfaces
with $m_{\text{pol}},n_{\text{tor}}=10$, which results in a smaller 661
dimensional optimization problem for each surface. See section III.4 for a
detailed description of the design decisions made for the code to handle the
complexity resulting from the bilevel optimization problem solved here.

The algorithm that we have followed in these experiments is:

  1. 1.

Obtain an initial coil set from a FOCUS-likeZhu  _et al._ (2018); Wechsung
_et al._ (2022) or near-axis expansion (NAE) optimizationGiuliani  _et al._
(2022b).

  2. 2.

Optimize the coils from step (1) for QA and nested flux surfaces using a
number of BoozerLS surfaces.

  3. 3.

Optimize the coils from step (2) using the BoozerExact surfacesGiuliani  _et
al._ (2022a) almost everywhere and keep BoozerLS surfaces only in the
neighborhood of troublesome low-order rationals.

As will be shown in the following examples, the magnetic field and physics
properties of the stellarator change drastically during step (2). For example,
islands and chaotic regions may appear and disappear from one iteration to the
next. Despite this, the robust BoozerLS surfaces can still be computed and
used for coil optimization. At the end of step (2), we have reached a
configuration with nested flux surfaces and much improved quasi-symmetry. In
step (3), we switch to the BoozerExact formulation as we are in the
neighborhood of the optimizer and do not expect islands to appear anymore. For
the BoozerLS surfaces that remain, we still include the PDE residual penalty
term in the objective, while for the BoozerExact surfaces, it is not needed.

We analyze the physics properties of the stellarator designs before and after
optimization by examining plots of the non-QA ratio $\|B_{\text{non-
QA}}\|/\|B_{\text{QA}}\|$, Boozer residual $\|\mathbf{r}\|^{2}_{2}$,
rotational transform on surfaces in the toroidal volume. Each point on the
profile corresponds to a BoozerLS surface that is computed with a continuation
procedure. The surfaces on which quasi-symmetry is optimized are verified by
comparing their cross sections with Poincaré plots. The physics profiles are
generated by computing many BoozerLS surfaces through the volume and
evaluating the physics quantities on each surface. The physics plots are
accurate away from island chains (Figure 1) when computed using BoozerLS
surfaces. Even though the ratio $\|B_{\text{non-QA}}\|/\|B_{\text{QA}}\|$ can
be computed, it is not well defined in regions of chaos and island as it
relies on the assumption that nested flux surfaces exist in the underlying
magnetic field. Moreover, when the rotational transform is computed using
BoozerLS surfaces, we find that the $\iota$ profile varies smoothly through
island chains. This is in contrast to what is observed when the profile is
computed with field line tracing, where the rotational transform is constant
through island chainsLandreman, Medasani, and Zhu (2021).

###  IV.1 Island healing

Figure 5: Island healing: Poincaré plot ($\phi=0$) and cross sections where we
optimize for quasi-symmetry and nested flux surfaces. The red cross section
corresponds to the $\iota=2/5$ surface and the outermost blue surface has
aspect ratio approximately 6. In the initial configuration, we use regularized
BoozerLS surfaces for robustness. After 1,000 iterations of the outer coil
optimization, we switch to BoozerExact surfaces (blue) but keep the
regularized BoozerLS surface at the low order rational $\iota=2/5$ (red).
Figure 6: Island healing: physics quantities before and after island healing.
The vertical red line corresponds to where the profile passes through
$\iota=2/5$ and the vertical blue lines correspond to the other surfaces on
which quasi-symmetry and nested flux surfaces were optimized. These figures
were generated using least squares surfaces, with
$m_{\text{pol}},n_{\text{tor}}=15$. We find that even though the island width
has been greatly reduced, the extremum in the the non-QA ratio has not
completely disappeared. The rotational transform profile varies smoothly
through island chains when computed using BoozerLS surfaces, which is in
contrast to what is observed when it is computed with field line
tracingLandreman, Medasani, and Zhu (2021).

We begin with an initial equilibrium field presented in (Ref. Landreman,
Medasani, and Zhu, 2021). This equilibrium was obtained by completing a stage
I optimization for quasi-symmetry, with a target aspect ratio of 6. Quasi-
symmetry was optimized on a single surface at $0.5$ normalized toroidal flux,
and the rotational transform was fixed to $0.39$ on the magnetic axis and to
$0.42$ at the boundary. That forces the rotational transform profile to pass
through $2/5$, a low-order rational, resulting in an island chain.

We found a coil set for this equilibrium by solving a FOCUS-like stage II
optimization problem (Refs. Zhu  _et al._ , 2018; Wechsung  _et al._ , 2022)
to obtain an initial set of coils to reproduce the stage I equilibrium
described in the previous paragraph. The stage II optimization was launched 16
times, each with slightly perturbed initial guesses. Once optimality was
reached, the best performing coil set was chosen to initialize optimization
described next. Poincaré plots reveal a large island chain at $\iota=2/5$
(figure 5). The non-QA ratio and BoozerLS residual through the plasma volume
in this initial configuration are shown in Figure 6, and we observe an
extremum in both curves as the rotational transform passes through
$\iota=2/5$. This is expected as we are attempting to fit nested magnetic
surfaces through a region with a significant island chain. After 1,000
iterations of the outer coil optimization, the extremum is no longer visible
and the slope of the rotational transform profile reverses.

The presence of the BoozerLS residual term in the objective results in a
trade-off with the non-QA penalty term. Switching to the faster BoozerExact
formulation everywhere except the low-order rational surface allows the
optimizer to further improve the quasi-symmetry in the volume. After the
optimization, both the non-QA penalty and Boozer residual were further
reduced. We observe that the optimization successfully reduced the width of
the island chain. The rotational transform profile still passes through the
low order rational number $2/5$ (figure 6), and there still is a local
extremum in the non-QA ratio plot, but its magnitude and width have been
substantially reduced.

###  IV.2 Chaos healing

Figure 7: Chaos healing: Poincaré plot ($\phi=0$) and cross sections where we
optimize for quasi-symmetry and nested flux surfaces. The innermost and
outermost red cross sections correspond respectively to the $\iota=2/5$
surface and the additional aspect ratio 4 surface.  Figure 8: Chaos healing:
physics quantities before and after optimization. The greyed region of the
plot corresponds to the volume between the outermost blue and red cross
section in Figure 7. Since this is a region with chaotic field lines and
without nested flux surfaces, the physics quantities are not well defined. The
vertical red line corresponds to where the profile passes through $\iota=2/5$
and the vertical blue lines correspond to the surfaces on which we optimize
for quasi-symmetry and nested flux surfaces. For comparison, we also include
the results from the previous chaos healing example, where we only optimized
on a smaller volume with aspect ratio 6.

The magnetic field generated by the initial coil set in section IV.1 presents
nested flux surfaces, a significant island chain, and chaotic field lines in
lower aspect ratio regions. The goal of this example is to simultaneously heal
islands as well as increase the volume in which nested flux surfaces are
present. We do this by adding a ninth surface with aspect ratio approximately
4 to the optimization problem.

A Poincaré plot of the initial and optimized configurations, along with cross
sections of the BoozerLS surfaces, is shown in Figure 7. A discussion of the
coils can be found in section V. We also stress the fact that in the initial
configuration, the outermost surface passes through a chaotic regions of the
magnetic field. We are able to compute such a BoozerLS surface thanks to the
least squares framework adopted in this work. Computing a BoozerExact surface
through the same regions of the field would be difficult since clearly nested
surfaces do not exist.

In Figure 8, the physics properties of the initial and optimized stellarators
are presented. Note that the non-QA ratio and rotational transform are not
well defined for the outermost surface in the original configuration, which we
indicate by the fading line. After initial coil optimization using the
BoozerLS surface description, reasonable magnetic surfaces are available and
thus we continue the coil optimization with the BoozerExact surface
formulation. The resulting design has nested flux surfaces as can be seen
visually from the Poincaré plot (bottom of Figure 7), and reasonable non-QA
ratio also out to the outermost surface (Figure 8). Comparing with the
previous example where we optimized on a smaller volume, we find that
requesting precise QA on lower aspect ratio surfaces causes the quality of QA
to degrade slightly at the core. This is unsurprising as the smaller the
region on which quasi-symmetry is requested, the easier it is to find a
magnetic field with that property Landreman, Medasani, and Zhu (2021);
Giuliani  _et al._ (2022b). In the next example, we present a configuration
optimized for QA only on the magnetic axis, which only presents extremely good
quasi-symmetry at the core, and whose quasi-symmetry quickly degrades moving
away from the axis.

###  IV.3 Cold start direct coil optimization

In the previous two examples, we followed the classical two-stage optimization
procedure to obtain coils, which we then improved using the optimization
problem (10). In this final example we show that thanks to the BoozerLS
formulation, this two stage procedure is unnecessary and we can design
stellarators with nested surfaces and precise QA properties starting from flat
coils.

The computation of Boozer surfaces in a magnetic field is more straightforward
when a magnetic axis with some nested flux surfaces are present. Moreover, the
solution to (1) is not unique when $\iota=0$. To avoid these issues, we start
with an initial configuration obtained from a near axis expansion (NAE)
optimizationGiuliani  _et al._ (2022b) beginning from equispaced, flat coils
with zero current (Figure 9).

A)B)C) Figure 9: Cold start direct coil optimization: Shown in (A) are the
flat coils with zero current used as initialization. The coils obtained with
the near-axis expansion optimization are shown in (B). These coils are used as
initialization for the BoozerLS and BoozerExact optimization, which result in
the coils shown in (C). The optimized coils from a previous work (Ref.
Giuliani  _et al._ , 2022a) obtained with a different local optimization
method are visually indistinguishable from the coils in (C). The coil currents
of the two designs are remarkably close to each other as well (Table 1).

In this initial NAE optimization, the geometry of the magnetic axis and
electromagnetic coils are optimized for QA on the magnetic axis. The coil
design requirements we impose are similar to those in Section IV.1, except
that we constrain each coil to have the same length of
$L=$4.5\text{\,}\mathrm{m}$$ rather than constraining the total coil length to
be less than $L_{\max}=$18\text{\,}\mathrm{m}$$. We found that when doing the
latter, the coils tend to have lengths that are disparate from one another as
only physics properties on the magnetic axis (as opposed to in the volume) are
targeted. The rotational transform on axis is constrained to $\iota=0.42$ and
the mean radial position of the axis is constrained to be 1.

Figure 10: Cold start direct coil optimization: Poincaré plot ($\phi=0$) in
the configuration obtained from the near axis expansion and after BoozerExact
optimization. The cross sections (blue curves) in the initial and final
configurations correspond to the surfaces used in the BoozerLS and BoozerExact
phases of the optimization. Despite the presence of island chains in the
initial configuration, the BoozerLS optimization successfully replaces them
with nested flux surfaces with precise quasi-symmetry.  Figure 11: Cold start
direct coil optimization: physics quantities after optimization from cold
start.

The NAE optimization is robust, and from a cold start finds a configuration
with precise quasi-symmetry on axis. However, the quality of the quasi-
symmetry quickly degrades as one moves away from the axis (Figure 11), as also
observed in Ref. Giuliani  _et al._ , 2022b. The rotational transform on axis
in the optimized configuration decreases from 0.42 and crosses the low order
rational $\iota=2/5$, which results in an island chain (Figure 10).

Using the coil set from the NAE optimization, we use eight regularized
BoozerLS surfaces, which pass directly through the island chain. After 1,000
iterations of BFGS, the islands have been healed and the rotational transform
no longer crosses a low order rational. Thus, we no longer need to use
BoozerLS surfaces and switch all the surfaces to BoozerExact. Since now there
no longer is a trade-off between the non-QA penalty and the Boozer residual in
the objective, one can more effectively improve the quasi-symmetry in the
volume. It is also notable how different the magnetic fields from the NAE and
BoozerExact optimizations are: the Poincaré plots and rotational transform
profiles drastically change.

The initial and final coils are shown in Figure 9. To our surprise, the
optimized coils found in this work visually overlap with the ones computed in
Ref. Giuliani  _et al._ , 2022a which correspond to the precise QA
configuration from Ref. Landreman and Paul, 2022. In table 1, the coil
currents found in this work and those in Ref. Giuliani  _et al._ , 2022a are
provided, where the current in the first coil is normalized to 1. We observe
that the currents in the two designs are remarkably close to one another, even
though the optimization algorithm starts with completely different initial
coils.

| Coil 1 | Coil 2 | Coil 3 | Coil 4
---|---|---|---|---
this work | 1.0000 | 1.06831 | 1.4101 | 2.2341
Ref. Giuliani  _et al._ ,2022a | 1.0000 | 1.06370 | 1.4037 | 2.2508
Table 1: Coil currents in the final designs of Ref. Giuliani  _et al._ , 2022a
and this work, normalized so that the first coil’s current is 1. We observe
that they are remarkably close to one another (within 1%) even though the
optimization algorithms start with completely different initial coils.

Finding similar coils and magnetic field as in Ref. Giuliani  _et al._ , 2022a
from scratch indicates that the surface optimization algorithm presented here
is not restricted to stay in a neighborhood of the initial coil set and that
the problem we solve here might suffer less from multiple minima.

##  V Conclusions and Discussion

In this work, we have presented a general algorithm to robustly search for
magnetic configurations with nested flux surfaces and precise QA, and which
can be realized by coils. This algorithm is a robust extension of our previous
work in Ref. Giuliani  _et al._ , 2022a in the following sense: in Ref.
Giuliani  _et al._ , 2022a, we only considered improvements of already
optimized stellarator magnetic fields that had nested flux surfaces for a very
large fraction of the total volume; in contrast, the extended algorithm we
describe here performs well even when we use unoptimized magnetic fields as
initial guesses for our optimizations. The proposed approach first generates
an initial coil set from either a FOCUS-likeZhu  _et al._ (2018); Wechsung
_et al._ (2022) or NAE optimizationGiuliani  _et al._ (2022b). Next, we use
the least squares surfaces formulation presented in this work to perform a
physics optimization and reach an area of coil parameter space with nested
flux surfaces. This is done by fitting nested surfaces, even in regimes where
flux surfaces do not exist. This is facilitated by the addition of a surface
area regularization to avoid self-intersections. After this optimization,
which behaves robustly even in the presence of island or chaos, the
stellarator design is improved using the local optimization algorithm
described in Ref. Giuliani  _et al._ , 2022a. In sections IV.1 and IV.2, we
show how to heal both localized islands and widespread chaos. In section IV.3,
we illustrate how to handle cold starts and show that our algorithm can lead
to optimized magnetic configurations that are drastically different from the
initial one. Despite a completely different optimization algorithm and
starting from flat coils and zero currents, we arrive at a coil set that is
visually indistinguishable from the one obtained using the local algorithm in
Ref. Giuliani  _et al._ , 2022a. This hints that designing coils for QA on a
volume is burdened by fewer local minima than targeting QA only on
axisGiuliani  _et al._ (2022b), or in a FOCUS-like optimizationZhu  _et al._
(2018); Wechsung  _et al._ (2022). This may also be due to the relative
restrictive maximal coil length value $L_{\text{max}}$ we have chosen; a
larger value might result in more local minima. We plan on running more
experiments to study this. In Table 2, we provide a comparison of the
geometric properties of the stellarator designs presented in this work along
with the geometric properties of their coils. Consider first the island and
chaos healing coil sets where we optimized up to surfaces with aspect ratio 6
and 4, respectively. When including surfaces with lower aspect ratio in the
optimization, the coil-to-surface separation decreases, and coil complexity
(maximum curvature, mean-squared curvature) reduces. At the expense of less
precise quasi-symmetry, we obtain a design with lower coil complexity and
nested flux surfaces with lower aspect ratio. In all examples, the coil-coil
separation inequality constraint was not active due to the chosen value of
$L_{\max}=$18\text{\,}\mathrm{m}$$. For longer coils, one can expect it to
become activeGiuliani  _et al._ (2022a).

Coil set |  Aspect ratio |  Coil lengths |  Maximum
curvatures |  Mean-squared
curvatures |  Coil-coil
separation |  Coil-surf.
separation
---|---|---|---|---|---|---
Island healing |  5.99 |  $4.08,5.37,4.37,4.17$ |  $4.21,3.96,4.12,3.97$ |  $5.00,4.08,5.00,5.00$ |  $0.108$ |  $0.254$
Chaos healing |  4.09 |  $3.93,5.05,4.85,4.17$ |  $3.76,3.37,3.61,3.64$ |  $4.66,3.12,3.80,5.00$ |  $0.126$ |  $0.195$
Cold start |  5.99 |  $5.37,3.94,4.10,4.60$ |  $4.28,3.64,4.17,4.76$ |  $4.27,5.00,5.00,5.00$ |  $0.108$ |  $0.240$
Table 2: Comparison of geometric properties of the stellarators and their
coils. Coil-to-surface distance is computed with respect to the outermost
surface on which quasi-symmetry is optimized.

This work opens the door to a more wider search for stellarator coils that
produce nested flux surfaces with precise quasi-symmetry. As a preliminary
study, we use the same setup as in the cold start example (section IV.3),
except that we consider various target average rotational transforms in the
plasma volume, namely $\overline{\iota}=0.1,0.2,\ldots,0.6$. Using the same
initial coil set from a NAE optimization with $\iota=0.42$ on axis, our robust
coil optimization procedure successfully finds coil sets with comparable
quasi-symmetry for a wide range of rotational transform values (Figure 12),
though larger rotational transforms corresponded to slightly worse QA.
Moreover, for the range of rotational transforms considered, we observe that
magnetic shear generally decreases with increasing $\overline{\iota}$. The
scan of $\overline{\iota}$ completed here used a fixed length
$L_{\max}=$18\text{\,}\mathrm{m}$$ for the sum of the 4 independent coils. It
is known that the total coil length is a strong regularizer for the outer
optimization problem, and improved quasi-symmetry can be obtained by
increasing $L_{\max}$, as shown in Ref. Giuliani  _et al._ , 2022a. In this
preliminary investigation, we always used the same initial coil set obtained
from a single NAE optimization to illustrate our method’s robustness and the
ability to find coils for a variety of rotational transforms. Alternatively,
one might recompute the initial coil set from a NAE optimization for each
$\overline{\iota}$ used in the parameter scan.

Figure 12: Quality of quasi-symmetry for various average rotational transforms
$\overline{\iota}$. Shown of the top are the non-QA ratios. We are able to
find stellarators with comparable quasi-symmetry for a wide range of average
rotational transforms $\overline{\iota}=0.1,0.2,\ldots,0.6$. Shown on the
bottom is the deviation from the average rotational transforms. All runs use
the shortest total allowable coil lengths $L_{\max}=$18\text{\,}\mathrm{m}$$
from the study in Ref. Giuliani  _et al._ , 2022a. For longer coils, one can
expect improved quasi-symmetry, but this question is not explored in this
article.

Using the tools presented in this work, we plan on completing a more detailed
scan of target physics and coil design values: average rotational transform
$\overline{\iota}$, $n_{\text{fp}}$, number of coils, maximum coil length
$L_{\text{max}}$, aspect ratio, etc. There are no fundamental restrictions
that prevent extending this approach to other flavors of quasi-symmetry such
as quasihelical or quasipoloidal symmetry. We focused here on vacuum-field
equilibria as they are an important aspect of stellarator optimization, and
can be used to initialize other optimization studies without the curl-free
assumption. An important extension of this work is to apply these algorithms
to direct coil optimization for finite $\beta$ magnetic equilibria. Doing this
requires the ability to compute the total magnetic field, along with its first
and second spatial derivatives at a given point in space. Although not
investigated here, it would also be interesting to develop an adaptive
algorithm that modifies the number of collocation points and surface
parameters on surfaces depending on the complexity of the surface.

## Code availability and optimized configurations

Many of the tools described here are available in SIMSOPTLandreman  _et al._
(2021), which is a suite of stellarator design, optimization, and analysis
utilities:

https://github.com/hiddenSymmetries/simsopt

## Acknowledgements

The authors would like to thank the SIMSOPT development team. This work was
supported by a grant from the Simons Foundation (560651). AG is partially
supported by an NSERC (Natural Sciences and Engineering Research Council of
Canada) postdoctoral fellowship. In addition, AC and FW are supported by the
United States National Science Foundation under grant No. PHY-1820852, and AC
is supported by the United States Department of Energy, Office of Fusion
Energy Sciences, under grant No. DE-FG02- 86ER53223. This work was supported
in part through the NYU IT High Performance Computing resources, services, and
staff expertise. The authors also gratefully acknowledge helpful discussions
with Elizabeth Paul.

## References

  * Giuliani  _et al._ (2022a) A. Giuliani, F. Wechsung, G. Stadler, A. Cerfon, and M. Landreman, “Direct computation of magnetic surfaces in Boozer coordinates and coil optimization for quasisymmetry,” Journal of Plasma Physics 88, 905880401 (2022a).
  * Helander (2014) P. Helander, “Theory of plasma confinement in non-axisymmetric magnetic fields,” Reports on Progress in Physics 77, 087001 (2014).
  * Landreman and Paul (2022) M. Landreman and E. Paul, “Magnetic fields with precise quasisymmetry for plasma confinement,” Phys. Rev. Lett. 128, 035001 (2022).
  * Landreman  _et al._ (2021) M. Landreman, B. Medasani, F. Wechsung, A. Giuliani, R. Jorge, and C. Zhu, “SIMSOPT: a flexible framework for stellarator optimization,” Journal of Open Source Software 6, 3525 (2021).
  * Hudson  _et al._ (2002) S. R. Hudson, A. Reiman, D. Strickler, A. Brooks, D. A. Monticello, and S. P. Hirshman, “Free-boundary full-pressure island healing in stellarator equilibria: coil-healing,” Plasma Physics and Controlled Fusion 44, 1377–1382 (2002).
  * Hudson  _et al._ (2003) S. Hudson, D. Monticello, A. Reiman, D. Strickler, S. Hirshman, L.-P. Ku, E. Lazarus, A. Brooks, M. Zarnstorff, A. Boozer, G.-Y. Fu, and G. Neilson, “Constructing integrable high-pressure full-current free-boundary stellarator magnetohydrodynamic equilibrium solutions,” Nuclear Fusion 43, 1040–1046 (2003).
  * Zhu  _et al._ (2019) C. Zhu, D. A. Gates, S. R. Hudson, H. Liu, Y. Xu, A. Shimizu, and S. Okamura, “Identification of important error fields in stellarators using the Hessian matrix method,” Nuclear Fusion 59, 126007 (2019).
  * Landreman, Medasani, and Zhu (2021) M. Landreman, B. Medasani, and C. Zhu, “Stellarator optimization for good magnetic surfaces at the same time as quasisymmetry,” Physics of Plasmas 28, 092505 (2021).
  * Baillod  _et al._ (2022) A. Baillod, J. Loizu, J. Graves, and M. Landreman, “Stellarator optimization for nested magnetic surfaces at finite $\beta$ and toroidal current,” Physics of Plasmas 29, 042505 (2022).
  * Cary (1982) J. R. Cary, “Vacuum magnetic fields with dense flux surfaces,” Physical Review Letters 49, 276 (1982).
  * Dommaschk (1982) W. Dommaschk, “Finite field harmonics for stellarators with improved aspect ratio,” Zeitschrift für Naturforschung A 37, 866–875 (1982).
  * Cary (1984) J. R. Cary, “Construction of three-dimensional vacuum magnetic fields with dense nested flux surfaces,” The Physics of fluids 27, 119–128 (1984).
  * Hanson and Cary (1984) J. D. Hanson and J. R. Cary, “Elimination of stochasticity in stellarators,” The Physics of Fluids 27, 767–769 (1984), https://aip.scitation.org/doi/pdf/10.1063/1.864692 .
  * Cary and Hanson (1986) J. R. Cary and J. D. Hanson, “Stochasticity reduction,” The Physics of Fluids 29, 2464–2473 (1986), https://aip.scitation.org/doi/pdf/10.1063/1.865539 .
  * Lee  _et al._ (2022) B. F. Lee, E. J. Paul, G. Stadler, and M. Landreman, “Stellarator coil optimization supporting multiple magnetic configurations,”  (2022).
  * Dudt  _et al._ (2022) D. Dudt, R. Conlin, D. Panici, and E. Kolemen, “The DESC stellarator code suite part iii: Quasi-symmetry optimization,”  (2022).
  * Boozer (2019) A. H. Boozer, “Curl-free magnetic fields for stellarator optimization,” Physics of Plasmas 26, 102504 (2019).
  * Dewar, Hudson, and Price (1994) R. Dewar, S. Hudson, and P. Price, “Almost invariant manifolds for divergence-free fields,” Physics Letters A 194, 49–56 (1994).
  * Giuliani  _et al._ (2022b) A. Giuliani, F. Wechsung, A. Cerfon, G. Stadler, and M. Landreman, “Single-stage gradient-based stellarator coil design: Optimization for near-axis quasi-symmetry,” Journal of Computational Physics 459, 111147 (2022b).
  * Landreman and Boozer (2016) M. Landreman and A. H. Boozer, “Efficient magnetic fields for supporting toroidal plasmas,” Physics of Plasmas 23, 032506 (2016).
  * Zhu  _et al._ (2018) C. Zhu, S. R. Hudson, Y. Song, and Y. Wan, “Designing stellarator coils by a modified Newton method using FOCUS,” Plasma Physics and Controlled Fusion 60, 065008 (2018).
  * Wechsung  _et al._ (2022) F. Wechsung, M. Landreman, A. Giuliani, A. Cerfon, and G. Stadler, “Precise stellarator quasi-symmetry can be achieved with electromagnetic coils,” Proceedings of the National Academy of Sciences 119, e2202084119 (2022).