# The Factory Must Grow: Automation in Factorio

Kenneth N. Reid  0001-8654-2430 Michigan State UniversityEast LansingMI48823
<EMAIL_ADDRESS>,  Iliya Miralavy  Michigan State UniversityEast
LansingMI48823 ,  Stephen Kelly  Michigan State UniversityEast LansingMI48823
,  Wolfgang Banzhaf  0002-6382-3245 Michigan State UniversityEast
LansingMI48823 and  Cedric Gondro  Michigan State UniversityEast
LansingMI48823

(2021)

###### Abstract.

Efficient optimization of resources is paramount to success in many problems
faced today. In the field of operational research the efficient scheduling of
employees; packing of vans; routing of vehicles; logistics of airlines and
transport of materials can be the difference between emission reduction or
excess, profits or losses and feasibility or unworkable solutions. The video
game Factorio, by Wube Software, has a myriad of problems which are analogous
to such real-world problems, and is a useful simulator for developing
solutions for these problems. In this paper we define the logistic transport
belt problem and define mathematical integer programming model of it. We
developed an interface to allow optimizers in any programming language to
interact with Factorio, and we provide an initial benchmark of logistic
transport belt problems. We present results for Simulated Annealing, quick
Genetic Programming and Evolutionary Reinforcement Learning, three different
meta-heuristic techniques to optimize this novel problem.

Factorio, operational research, video games, optimization

††copyright: rightsretained††doi: 10.1145/nnnnnnn.nnnnnnn††isbn: 978-x-xxxx-
xxxx-x/YY/MM††price: 15.00††ccs: Theory of computation Convergence and
learning in games††journalyear: 2021††copyright: acmlicensed††conference: ;
15.00††price: ;††doi: 10.1145/3321707.3321769††isbn: 978-1-4503-6111-8/19/07

##  1\. Introduction

Factorio is a video game produced by Wube software111https://factorio.com/,
which was fully released in August 2020. The premise of Factorio is an
engineer who crash-lands into an alien planet, and must use their engineering
expertise to accumulate resources, automate production and launch a rocket to
escape the planet, all while avoiding the aliens defending their world from
the pollution created by the engineers machinery.

Games represent a wide variety of challenging and interesting problems
designed to engage human intelligence. As such, games provide a rich test bed
for Artificial Intelligence. Yannakakis and Togelius (yannakakis2018, )
structure the application of AI in games into three categories: 1) AI for game
playing; 2) AI for content generation; and 3) AI for player modelling. Of
particular interest to this work are the first two categories on which we
shall focus from here on. Of particular interest to this work are the first
two categories on which we shall focus from here on.

In AI for Game Playing, artificial agents learn to interact with a game,
observing the state of play and taking actions in an attempt to maximize one
or more objectives. In general the objective of an artificial game-player may
be to win the game (e.g. maximize game score or length of play) or to define
behavior of a Non-Player Character (NPC) in the game such that the game
experience is more engaging for a human player. Game playing is the oldest and
likely the most common application of AI in games, and a wide array of
literature exists on the subject. Some recent and well-known examples include
artificial soccer players in RoboCup (Asada2019, ), mastering the board game
Go with Artificial Neural Networks (alphazero, ), and playing classic Atari
games at a level comparable to human players (mnih15, ). Artificial game
playing agents have many practical applications, from discovering bugs in
computer games or other interactive software (chrabaszcz2018, ), to autonomous
robotics (Dickinson2007, ).

Procedural Content Generation (PCG) refers to the use of AI to create game
content automatically, with minimal human input. Game content may include
environments, difficulty levels, rules, or characters. There are many
compelling reasons to employ PCG. When game design is automated, one could
imagine building a game that never ends, with PCG continuously creating new
content. PCG could also support games tailored to individual players, with the
goal of improving the entertainment value and/or educational value of a game
(Cruz2017, ). In short, PCG extends the creativity of human game designers and
can therefore be used to study computational creativity more generally, see,
e.g., (Phillips2016, ).

Factorio provides a breadth of opportunities for AI Game Playing and PCG that
directly transfer to real-world problem solving scenarios. To simplify the
complex mass of problems facing the player, herein mentioned problems include
the maximization of the players research progress; the efficient usage of
resources, space, the minimization of pollution, and more. In particular, the
link between Operational Research (OR) and environmental issues is of a
potential use of Factorio (bloemhof1995interactions, ). Factorio is
essentially an open-ended collection of operational research and optimization
problems.

Much like an ecology, Factorio can be thought of as a large set of varied
interlinking problems and solutions, where the key to success is often
balancing between solving these problems. For example, focusing solely on
production will likely increase the pollution output, causing a game failure.

In this paper, we provide an overview of some such problems and provide an
interface by which users can connect external problem solvers to Factorio. As
an example of this interface we tackle a small but common set of problems: the
logistic belt problems. Belts transport raw resources and manufactured goods
between collection points, machinery and storage. With hundreds of items in
the game which can be delivered this way, often requiring multiple belts
dedicated to each item, this is both a complex and common problem facing
players. Additionally, there are different types of belts, each with different
potential speeds. The purpose of these belts is to deliver goods, but also to
meet specific criteria: the correct ratio of items output by the solution;
materials reaching the correct end destinations; belts avoiding obstacles; and
minimizing the number of belts used in a solution.

In this paper we solve several logistic belt problems of varying grid sizes
and number of obstacles, by employing various metaheuristics through the
interface to test the generated solutions within Factorio itself. Solving
problems which are analogies for real-world problems is beneficial …

##  2\. Related Work

There are several works in other fields which mention Factorio: virtual
economy analysis (mizutani2019toward, ), video game learning and community
building (kunze2019video, ), and a consideration of the capitalistic and
production-obsessed players of Factorio in (leblanc2018only, ). In a poster by
Duhan et al (duhan2019factory, ) the concept of creating a reinforcement
learning based engine to play Factorio is described, with the potential for
creation of transferable knowledge to real-world scenarios, though no further
work has yet been published.

In Factorio logistic belts are used to solve problems of transporting
materials. A clear similarity between this contribution and real-world
problems is the use of conveyor belts, their purpose, functionality and
optimization. Conveyor belts, in the real world, have been referred to as
productive and a particularly economical transport mechanism with high
performance compared to alternatives and a small carbon footprint
(grinvcova2014experimental, ). A recent survey paper details the importance of
Automated Sortation Conveyor (ASS) (boysen2019automated, ) \- highlighting
their efficiency, time-savings and the reduction in employee costs, with more
details on the importance of logistics for transporting items via belts. The
overland conveyor is often used in transporting raw materials; such as ore
before processing. Much of the literature on the optimization of such devices
is with regard to design of the machine itself, and the energy consumption of
fully-loaded and empty conveyors, rather than planning logistical routes
(gerard2009optimisation, ) (carr2020reducing, ). Particle Swarm Optimization
(PSO) was also used in a similar optimization of screw conveyors, a type of
conveyors used to transport granular bulk goods, but also considered design of
the machinery rather than logistical placement and application
(cai2011optimisation, ).

Routing problems are relevant to this paper as the problem being solved also
requires a connection between two points on a grid. A masters thesis by one of
the Wube Software developers, Ondřej Majerech (majerech2017solving, ), which
considers the efficient path planning of numerous agents in a dynamic world,
is one such problem. While our routing problem is not constrained by problems
such as time, the agents in the world require efficient algorithms to swiftly
plan routes, to prevent players from encountering enemies which are static. In
the real-world, there are multiple problems which require optimization to find
feasible or better than feasible solutions. Vehicle routing problems, first
defined by Dantzig et al (dantzig1959truck, ), is similar to other complex
problems such as the Travelling Salesman Problem (TSP) in that an agent is
placed at a starting point in a network and, following some constraints, must
reach a certain goal. Constraints such as backtracking restrictions or maximum
visits to each location, or even additional stops, can be applied. In the
survey paper by Pillac et al (pillac2013review, ) a large number of sub-
domains are presented, including the transportation of people, goods and
services, where the need for optimization is clear. Additional real world
problems for routing can be found in the expanding domain of network routing.
Such research areas include Device to Device (D2D) network routing
(shaikh2018routing, ), unmanned aerial vehicle networks (arafat2019routing, )
and disaster relief networking (jahir2019routing, ). Each of these research
areas have differing constraints, with for example, disaster relief having
requirements for redundancy and flexibility that other network optimizations
may consider less important.

Similar to routing problems is circuit planning, due to the need to prevent
overlapping signals, which is analogous to preventing transport belts from
passing materials from one to the other. An interesting outcome from a
landmark paper in circuit design optimization using metaheuristics required a
user-supplied cost function, placing the onus of constraint enforcement and
relaxation on the user (gielen1990analog, ). This is similar to how our
constraints are provided; before analysis has begun on feasibility. A more
recent problem of circuit design and efficiency used PSO, an effective
algorithm for non-linear exploration of a problem space, to discover efficient
pathways for signals (fakhfakh2010analog, ). Analog block design is
qualitatively evaluated from a perspective of maximizing performance, rather
than minimizing cost, though the two may be causally related.

##  3\. Problem Description

There are a variety of problems in Factorio which can be automated and
optimized through the use of metaheuristics. For this initial work, we study a
set of problems from the perspective of direct optimization and machine
learning. These methods differ with respect to their goal. Direct optimization
searches in the space of solutions directly for a single point which minimizes
a problem-specific cost function. Machine learning, on the other hand,
optimizes a prediction machine which accepts a representation of the problem
as input and extrapolates, or predicts, the optimal solution point relative to
the input state. As such, the goal of direct optimization is to discover a
single optimal solution relative to a finite training data-set (Section 5.2),
while the goal of machine learning is to discover a solution that generalizes
to unseen data (Sections 5.3 and 5.4). Put another way, direct optimization
builds solutions, while machine learning builds problem solvers. The later
approach potentially incurs a significantly greater computational cost because
a larger number of training examples, and therefore model evaluations, are
required to build generalized behaviors.

In this paper, we focus on the logistics belt placement problem. In order to
best represent the various aspects of belt placement, we define six problem
instances with varying attributes. The problem matrices provide empty cells in
sizes of $3\times 3$, $6\times 6$ and $12\times 12$. Outside of these empty
matrices are walls, two cells thick, an inserter (for input) with a chest, and
a receiver (for output) with a chest. For a visual representation, see the in-
game representation in Fig 1, where the input inserter and chest are at the
top right of the image, and the output receiver and chest are at the bottom
left. Underground belts function by requiring an input and an output belt to
be in the same row or column, face the same direction, and have a maximum of 4
squares between the two items. Without this criteria being met, underground
belts become nothing more than obstacles.

Figure 1. 12 x 12 matrix with obstacles

As well as the three aforementioned empty grids, we tackle three other problem
matrices which are identical in size but have obstacles positioned where
logistic items cannot be placed. This both simplifies and complicates the
problem, as there are fewer possible solution combinations, with some of them
no longer able to transport materials from source to destination.

For accessible manipulation by external optimizers, the solutions and problems
are encoded as integer matrices. The interface, described in Section 5,
translates the integers (representations of which are shown in Fig 2) and
provides commands to Factorio to create the problems and solutions in-game for
experimentation. Integers less than zero are not utilized by solutions, but
instead are used to represent the problem.

Figure 2. Integer Encodings

The quality of a solution is measured through its fitness value. A fitness of
0 is assigned if a solution is deemed infeasible (having failed one or more
hard constraints (HCs), or a value between 0 and 1 otherwise. The hard
constraints can be listed as follows:

  * •

HC1: A maximum of one logistic object can be placed in a single position at a
time or one obstacle.

  * •

HC2: There must be at least one logistic item placed for a solution to be
deemed feasible.

If a solution is feasible, then each Soft Constraint (SC) can be evaluated to
return a single fitness value from the objective function. The soft
constraints can be listed as follows:

  * •

SC1: Maximize the number of item output by the solution.

  * •

SC2: Maximize the number of item input into the solution.

##  4\. Mathematical Model

This section of the paper concisely defines the problem in a mathematical
model, though each algorithm in this paper has a unique implementation. Table
1 provides the meaning for the mathematical objects in use.

Parameters |  Description
---|---
$M$ |  A parameterized matrix of two vectors, $M_{x}\in M$ containing all $x$ coordinates and $M_{y}\in M$ containing all $y$ coordinates. Coordinates (position) of an individual position are referred to as $p$, which is a vector with a single $x$ coordinate and $y$ coordinate.
$d$ |  Cardinal direction $d\in\\{0,1,2,3\\}$ which represent north, east, south and west, respectively.
$b_{pd}$ |  A belt which can receive 0 or 1 inputs, and 0 or 1 outputs.
$u_{pd}$ |  An underground belt which can receive 0 or 1 inputs and 0 or 2 outputs.
$L$ |  An array of two vectors containing all belts and underground belts.
|  $b\in L_{b}$, $u\in L_{u}$
$w_{k}$ |  Weightings for each SC.
|  $k\in\\{1,2\\}$.
|  $w_{k}\in\mathbb{R},~{}~{}k\in\\{1,2\\}$
|  $0\leq w_{k}\leq 1,$ $k\in\\{1,2\\}$
|  $w_{1}+w_{2}=1$
$\phi_{p}$ |  Obstacle at position $p$. Equal to 1 if there is an obstacle, 0 if not.
|  $\phi\in\Phi$.
$i$ |  A parameterised number of items to be input to the solution.
|  $i>0$, $i\in\mathbb{N}$
$o$ |  A parameterised number of items to be output by the solution.
|  $o>0$, $o\in\mathbb{N}$
Table 1. Parameters

Decision variables are the mechanism for directing the solution towards a
higher quality and feasible solution. Decision variable (1) is used for
placing a belt at any position $p$ in the provided grid $G$, of any type $c$
and of any direction $d$. A belt can either be placed, $1$, or not, $0$.
Similar decision variables are required for underground belts $u$, equation
(2).

(1) $\displaystyle X_{b_{cpd}}=0~{}\textrm{or}~{}1\textrm{ , }$
$\displaystyle\forall b\in L_{b},$ $\displaystyle\forall p_{x}\in
G_{x},\forall p_{y}\in G_{y},$ $\displaystyle\forall d\in\\{0,1,2,3\\}$

(2) $\displaystyle X_{u_{cpd}}=0~{}\textrm{or}~{}1\textrm{ , }$
$\displaystyle\forall u\in L_{u},$ $\displaystyle\forall p_{x}\in
G_{x},\forall p_{y}\in G_{y},$ $\displaystyle\forall d\in\\{0,1,2,3\\}$

The objective function is defined as maximizing the sum of each weighted SC.

(3) $\displaystyle\textrm{Max:}~{}~{}(SC1+SC2)$

Subject to:

(4) HC1: $\displaystyle\Bigg{[}\phi_{p}+\Bigg{(}{\sum}_{b\in
L_{b}}X_{b_{p}}\Bigg{)}+\Bigg{(}{\sum}_{u\in
L_{u}}X_{u_{p}}\Bigg{)}\Bigg{]}<2~{}~{}$ $\displaystyle\forall p_{x}\in
M_{x},\forall p_{y}\in M_{y}$ $\displaystyle\forall p1_{x}\in M_{x},\forall
p1_{y}\in M_{y}$ $\displaystyle\forall\phi\in\Phi$

(5) HC2: $\displaystyle\Bigg{[}\Bigg{(}{\sum}_{b\in
L_{b}}X_{b}\Bigg{)}+\Bigg{(}{\sum}_{u\in L_{u}}X_{u}\Bigg{)}\Bigg{]}>0~{}~{}$

(6) $\displaystyle\textrm{SC1:}~{}~{}o*w_{1}$

(7) $\displaystyle\textrm{SC2:}~{}~{}i*w_{2}$

##  5\. Implementation

The source code of Factorio cannot be modified directly; instead there are
three mechanisms where users can provide input to the game. These are 1)
through mouse and keyboard control, which is the standard use case for
Factorio players, 2) through commands delivered via the in-game console, 3)
via ‘mods’ (short for modifications). Direct usage of none of these mechanisms
are sufficient for an optimizer to interact with Factorio and evaluate
possible solutions to in-game problems. In section 5.1 we introduce Factorio
Optimizer Interface (FOI), an interface which uses Remote Console (RCON)
protocol and modding to connect to the Factorio in-game console and enables
optimizers, regardless of their programming language, to easily interact with
a running Factorio server.

Modding (the act of writing mods) is often supported by video game companies.
With Factorio, the company has provided a modding interface through the
scripting language of Lua. While restricted to prevent run-time errors and
modification of core elements of the Factorio software, Lua provides
structured access to modifying many elements of Factorio. For example, many
mods are written to add graphical elements, music or sound effects, additional
game-play features such as new vehicles or modifying how existing mechanics
function.

RCON (rconsource, ) is a TCP-based protocol which uses a request-response
communication method to enable server conductors to remotely administrate
servers. RCON establishes a connection by using a simple authentication method
which includes sending a password to a server and receiving a response. Upon
success, a connection is established and the conductor is able to communicate
with the server by sending request packets and receiving responses
(rconsource, ).

###  5.1. Factorio Optimizer Interface

We introduce FOI, an open-source cross-platform interface written in Python
that uses RCON protocol to enable problem solvers to connect and interact with
a running Factorio server console. FOI incorporates file system and command-
line interface, two programming language-independent methods, to communicate
with problem solvers.

Figure 3 illustrates an optimization system that uses FOI to interact with the
Factorio server. In each iteration, the optimizer prepares matrices
representing the problem and solution(s) (FOI is capable of evaluating
multiple solution matrices at the same time) and calls the interface program
giving it the path to the matrices as input. FOI then, uses integer encodings
(Figure 2) of the matrix cells to create corresponding Lua commands to place
entities on the game grid as it is explained in section 3. The RCON protocol
enables FOI to send these Lua commands to the game server to execute them and
to receive the server response. For the logistic belt placement problem,
server response will show the number of resources taken from the input chest
by the input inserter and the number of resources successfully transferred to
the output chest. FOI sends the server response back to the optimizer so the
fitness of the solution matrix can be evaluated using these information. Some
entities in Factorio require certain elements to become enabled. For instance,
an inserter needs electricity to transfer resources from/to transport belts.
Designing power grids is a cumbersome problem in Factorio. To avoid solving
that issue while building problem matrices, FOI uses the interface mod to
remove such constraints and minimizes the amount of entities needed to be
placed on the Factorio grid for a problem to be running in-game while
remaining visually informative.

Figure 3. Diagram of FOI

In this paper we present a number of metaheuristics which, through the modding
interface provided by Wube software and RCON, test solutions to the problem.
The selected algorithms represent a variety of heuristic types, including
reinforcement learning, genetic programming and classical metaheuristic
techniques. Due to the universal inputs allowed by the interface presented in
this paper, the metaheuristics are implemented in Java, Python and C++. As
described in Section 1, the complexity of the problem can be expanded to
include different speeds of logistic machinery, additional belt-types (such as
splitters which can be used to pass objects from one belt to many, many to
many, or many to one). As such, the selected algorithms which are not simple
path-finding algorithms in order to establish an initial benchmark for this
combinatorial optimization problem.

###  5.2. Parallel Simulated Annealing

Simulated Annealing (SA) is a classic search heuristic in the history of
metaheuristics. This metaheuristic was initially described in 1953 by
Metropolis et al (metropolis1953equation, ) as an adaptation from the
Metropolis-Hastings algorithm. This algorithm was selected partly due to it
its popularity, allowing this to be a useful insight into this problem and how
it can be solved. Additionally, we selected this algorithm due to its usage in
real-world problems, such as mobile network routing (saha2017seer, ), vehicle
routing (vincent2017simulated, ), systems-on-chip design flow
(gavrilov2018clustering, ).

We utilize a version of SA with simultaneous solution comparisons known as
Parallel Simulated Annealing (PSA). This is largely due to the inefficient
nature of a serial SA within the constraints of this problem, as well as due
to the proven ability of similar algorithms such as Differential Evolution to
solve larger search spaces, for example by Whalen et al (whalen2020evolving,
). In order to facilitate this, while still maintaining the core attributes of
SA, the algorithm works as can be seen in the Algorithm 1.222The variable
lettering in this section is unrelated to other sections of this paper.

The number of iterations the PSA algorithm will perform is calculated as
$\big{\lceil}\log{(\delta/t)}/\log{(cr)}\big{\rceil}$. The selection of a
neighbor in line 1 is conducted by modifying a bounded-random number of
entities, based on the current number of iterations compared to the total
number of iterations to be carried out. This allows the potential number of
changes to the solution integer array to be bounded highly at the beginning of
run-time. As the temperature $t^{\prime}$ reaches $k$ the number of changes to
select a neighbor will be much lower, with a lower bound of one.

Fitness is calculated in lines 1 and 1 by evaluating each solution in the
vectors, and enumerating inputs and outputs each solution allows.

Following traditional SA functionality, solutions are always accepted if
fitness is improved, and may be accepted if fitness worsens, depending on run-
time and potential fitness loss. The method this decision is described in the
function at line 1 $n$ times.

Data: $t$ Initial temperature

Data: $t^{\prime}$ Current temperature

Data: $\delta$ Delta

Data: $cr$ Cooling Rate

Data: $n$ Solutions per iteration

Data: $\vec{S}$ Previous solutions

Data: $\vec{S^{\prime}}$ Current solutions

Data: $\vec{S^{\prime\prime}}$ Temporary solution vector

1

2$\vec{S}\leftarrow stochasticallyGenerateNSolutions(n)$

3

4while  _$t^{\prime} >\delta$_ do

5 for (  _$i\leftarrow 0;\ i <n;\ i\leftarrow i+1$_ ) $\\{$

6 $\vec{S^{\prime}}_{i}\leftarrow selectNeighborOf(\vec{S}_{i}$)

7

8 $\\}$

9 $\vec{F}\leftarrow calculateFitness(\vec{S})$

10 $\vec{F^{\prime}}\leftarrow calculateFitness(\vec{S^{\prime}})$

11

12 for (  _$i\leftarrow 0;\ i <n;\ i\leftarrow i+1$_ ) $\\{$

13 if  _$\vec{F^{\prime}}_{i} >=\vec{F}_{i}$_ then

14 $\vec{S^{\prime\prime}}_{i}\leftarrow\vec{S^{\prime}}_{i}$

15

16 end if

17 else if
_$acceptanceProbability(t^{\prime},\vec{F^{\prime}}_{i},\vec{F}_{i})$_ then

18 $\vec{S^{\prime\prime}}_{i}\leftarrow\vec{S}_{i}$

19

20 end if

21

22 $\\}$

23 $\vec{S}\leftarrow\vec{S^{\prime\prime}}$

24 end while

25

26 Function  _acceptanceProbability(_$t^{\prime}$ ,$f^{\prime}$,$f$_)_:

27 $r\leftarrow randomDouble()$

28 $b\leftarrow 1-f$

29 $b^{\prime}\leftarrow 1-f^{\prime}$

30

31 if  _$r >e^{b^{\prime}-b}/t^{\prime}$_ then

32 return _true_

33

34 end if

35 return _false_

36

37

38

Algorithm 1 Parallel Simulated Annealing

###  5.3. qGenetic Programming

Genetic Programming (GP) is an evolutionary algorithm which evolves computer
programs to solve a problem that can be computationally represented. Initially
a population of individuals is created and their fitness is evaluated. The
best individual(s) are preserved and a portion are selected based, determined
by fitness, to produce the next generation of individuals using evolutionary
operators such as crossovers and mutations (Koza1992, ).

We present q Genetic Programming (qGP) for the first time in this paper to
solve the logistic belt problem (algorithm 2). This algorithm represents the
solutions with a linear genome consisting of a mixture of strings and integers
that can be mapped to operators (functions or non-terminals) and operands
(terminals). The goal is to interpret each genome to connect as many operators
as possible in the order of their appearance with suitable operands to form a
competitive model. For instance, a simple summation operator requires two
numeric operands to produce a meaningful output. The output of each function
then, can be considered as an operand for the next possible operators. qGP
neglects occurrence of an operator in case the algorithm fails to find fitting
operands for it. A sequence of connected operators form a solution or a model.

In order to interpret a genome qGP uses a mapping table to translate each base
value to an operator or an operand. There is no limit to the data types that
can be used as operands in qGP. In fact, it is designed so new data types for
operands and/or operators with unique functionalities can be easily added to
the optimizer. qGP creates a list of all the operands and operators in the
order they appear in the genome. Starting from the first operator, if the
algorithm can satisfy the operator’s operand requirements, then it will use
the operator on a set of applicable operands. The operator will then be output
to the end of the operands list. If the operator demands cannot be satisfied,
qGP saves the operator in a different list called idle operators. If a new
value is added to the operands list, idle operators become eligible for
consideration again to check if the new operand will enable them. In qGP,
unlike many GP representations, operators can have multiple outputs or inputs
which allows it to build complex model structures. Additionally, outputs of
operators can be used multiple times with a single occurrence in the genome.
This promotes modularity in qGP models while keeping the genome as compact as
possible. This is similar to gene duplication in nature but without any
overhead in the genome. After interpreting models into combinations of
operators and operands (in our case, sequence of calls to Connect() and
Place() functions), qGP saves the matrices generated from each individual to
give to FOI for fitness evaluation.

Table 2 lists the parameters and configurations used in this paper to
undertake the logistic belt problem using qGP. We use two operators named
Place and Connect in our function setPlace requires a matrix, a grid position,
and an object (a belt or underground belt) as an input. The operator then
places the object on the given position of the input matrix and returns the
enhanced matrix as an output. The Connect operator requires two grid positions
and a matrix as input to logistically connect the two positions using
transport belts. Connect can place multiple objects on the grid simultaneously
while Place only assigns one object on a grid at a time. We use Tournament
Selection (Blickle1995, ) of size five and elitism (Poli2008, ) of one
individual in forming the next generation of solutions.

Parameter |  Value
---|---
Population Size |  30
Initial/Max Genome Size |  40/100
Tournament Size |  5
Elitism |  1 individual
Operators |  Connect, Place
Operands |  4 belts and 4 underground belts for each direction, problem matrix, integers (limited to problem matrix size)
Table 2. qGP Parameters

Data: ${g}$ Current generation

Data: ${G}$ Maximum number of generations

Data: ${P}$ Population of solutions

1 ${P}\leftarrow initializePopulation()$

2 $best\leftarrow 0$

3 while  _$g <G$_ do

4 for (  _$i\leftarrow 0;\ i <size(P);\ i\leftarrow i+1$_ ) $\\{$

5 $interpretModel(P[i])$

6 $createMatrices(P[i])$

7 $evaluateFitness(P[i])$

8 if  _P[i].fitness > best_ then

9 $best\leftarrow P[i].fitness$

10

11 end if

12

13 $\\}$

14 $newPop\leftarrow[]$

15 while  _size(newPop) < size(P)-1_ do

16 $\vec{selected}\leftarrow tournamentSelection(P)$

17 $newPop\leftarrow newPop+crossover(\vec{selected})$

18

19 end while

20 $newPop\leftarrow mutate(newPop)$

21 $newPop\leftarrow newPop+best$

22

23 end while

Algorithm 2 qGenetic Programming

###  5.4. Evolutionary Reinforcement Learning

Reinforcement Learning (RL) is a type of machine learning that simulates the
way living organisms adapt through trial-and-error interaction with their
environment over time(sutton18, ). When modeling this interaction in a
computer simulation, time is divided into discrete steps. At each timestep,
the digital organism, or agent, observes its environment, takes an action that
changes the state of the environment, and receives a feedback signal that
describes the desirability of its current situation. The goal in RL is to
develop agents that map observations to actions such that the summed feedback,
or reward, over all timesteps is maximized. Most real-word decision-making
problems can be modeled in this way, including game playing. In this work,
Factorio provides a rich ecosystem of problem environments, and agents are
represented by computer programs. This Section describes how we can use GP in
combination with Reinforcement Learning to evolve programs that solve multiple
problems in Factorio by interacting with the game. In brief, a population of
programs will be stochastically initialized, evaluated in the Factorio
environment, selected for quality, and modified through an iterative process
inspired by the Darwinian principle of natural selection. This Evolutionary
Reinforcement Learning (ERL) process is summarized in Algorithm 3.

Interactions with a video game are naturally episodic. The initial state of
the environment is set by the game engine and the player interacts with the
game until a terminal state is reached (e.g. Game Over) or a time constraint
is exceeded. Human players typically need to play multiple episodes in order
to gain the breadth of experience required to become a skilled player. Being
naturally curious, they may try something new in each episode, leading to a
diversity of learning experiences. By contrast, our game-playing programs
behave deterministically: given the same initial conditions and a
deterministic environment, the programs will produce the same sequence of
actions in every episode. Thus, to ensure the programs are exposed to a
variety of scenarios, each program is evaluated in multiple episodes with
unique environmental conditions (See Algorithm 3, lines 3 to 3). In this work,
each episode represents a unique belt problem. Then $genProblem$ function in
Algorithm 3 randomly generates problem matrices of a given size, with or with
obstacles. The quality, or fitness, of a program is measured by its average
reward over all episodes it experiences during its lifetime. Thus, programs
that perform well in as many belt problems as possible will be favored by
selection operators. Programs are evolved using a generational Genetic
Algorithm (GA). In each generation, every program is evaluated in 5 unique
episodes, up to a lifetime maximum of 100 episodes per individual (that is, if
they survive for 20 generations). After all programs are evaluated, they are
sorted by fitness and the worst 50% are deleted from the population (Algorithm
3, line 3). New programs are created by uniformly sampling two parents from
the surviving individuals and applying crossover and mutation operators. This
processes is repeated until the population is back to $P_{size}$(Algorithm 3,
lines 3 to 3). Over time, this combination of GP and RL develops generalized
game-playing programs capable of solving multiple unique belt problems in
Factorio. The particular GP framework we employ for this study is Tangled
Program Graphs (TPG). No problem-specific functions were used. A more detailed
account of TPG in RL and video game playing can be found in (kelly18a, ).

####  5.4.1. RL Interface and Evaluation

The observation space in our formulation of Factorio belt problems takes the
form of an integer matrix (See Figures 2 and 1). However, the RL agents in
this work operate from a partial world-view. At any timestep, the agent is
situated in one cell of the problem matrix and its observation space,
$\vec{s}(t)$, consists of the 3x3 matrix of cells centered at the currently
occupied location, for a total of 9 integer state variables in $\vec{s}(t)$.
The agent is initialized in the center of the matrix at the beginning of each
episode. The location of inputs and outputs will never be visible from this
initial state, and the agent must therefore explore its surroundings and
possibly retrace its steps when deciding on a path for transport belts333The
TPG variant employed in this work includes a temporal mechanism (kelly2020, )
which is mostly likely useful in the partially-observable state space, but an
ablation study to confirm the significance of memory is left to future work..
In each timestep, the agent can choose to move one step in any cardinal
direction, or place one of 4 types of transport belt at the current location.
Thus, thus total action space is $\mathcal{A}={N,E,S,W,1,2,3,4}$,
corresponding to one-step movement in any direction or the placement of one
type of transport belt at the agent’s current location. Transport belts can be
placed in any cell to overwrite the previous cell contents, except where this
would overwrite a component of the initial problem (i.e. walls, chests, or
inserters). Naturally, the agent’s location is limited to the bounding box
defined by the current matrix size (3x3 or 6x6 in this work). Each episode
lasts for 20 timesteps. When this time constraint is reached, the final matrix
state is submitted to Factorio for evaluation.

Data: ${g}$ Current generation

Data: ${G}$ Total number of generations

Data: ${P}$ Population of programs

Data: $P_{size}$ Program population size

Data: ${M}$ Set of problem matrices

1 $P\leftarrow initializePopulation()$

2 while  _$g <G$_ do

3 $M\leftarrow generateProblems$

4 for (  _$i\leftarrow 0;\ i <size(P);\ i\leftarrow i+1$_ ) $\\{$

5 for (  _$j\leftarrow 0;\ j <size(M);\ j\leftarrow j+1$_ ) $\\{$

6 $evaluateFitnessRL(P[i],M[j])$

7

8 $\\}$

9

10 $\\}$

11 $P\leftarrow deleteWorst(P)$

12 while  _$size(P) <P_{size})$_ do

13 if  _$randomDouble() <crossoverProbability$_ then

14 $parent1\leftarrow sampleUniform(P)$

15 $parent2\leftarrow sampleUniform(P)$

16 $child\leftarrow crossover(parent1,parent2)$

17

18 end if

19 else

20 $child\leftarrow sampleUniform(P)$

21

22 end if

23 $P\leftarrow P+mutate(child)$

24

25 end while

26

27 end while

Algorithm 3 Evolutionary Reinforcement Learning

##  6\. Results

We utilize different metaheuristics presented in section 5 to investigate the
performance of each algorithm in this novel problem. Testing was conducted on
mixed hardware, but each heuristic was given approximately equal maximum
computational time by limiting number of iterations. While the authors
recognize that an iteration in one heuristic is not equivalent to an iteration
in a differing heuristic, the purpose is to investigate varying algorithm
types rather than compare similar heuristics and report the best.

(a) 3x3 Problem

(b) 6x6 Problem

(c) 12x12 Problem

Figure 4. Average Best-Fitnesses by Iterations. (y-axis of 6x6 and 12x12 plots
are log scale)

Figure 4 plots the average best fitness of the best solution discovered by
each algorithm over 167 iterations for each problem matrix. qGP provides the
highest scoring solutions in all cases but one. This is likely a result of the
powerful task-specific operators included in the qGP function set (See Connect
Section 5.3 Table 2). However, in the 6x6 problem the obstacles caused the
task-specific operators to fail to find a high quality solution. PSA
successfully finds solutions for the 3x3 problems but cannot find solutions
for the 6x6 nor the 12x12 problems. PSA has similar performances with and
without obstacles, whereas qGP shows distinctly different attributes when
solving obstacles compared to no obstacles for each of the three problems. In
the 3x3 problem (and the 6x6 and 12x12 for qGP), the initial steep climb in
fitness is due to solutions being created which allow transportation of
objects from the input without them being delivered to the output. Additional
success once the two points are connected can be due to more efficient
pathways allowing more objects to be delivered before allowed run-time
completes. Fitness in the solutions such as those produced by PSA in the 6x6
problem, for both with obstacles and without, is continuing to climb as the
distance covered by additional belts is increased, however they fail to reach
the output. Due to the high number of possible changes that can be made (9
integers to select from in each cell, is $9!^{12*12}$), SA fails to compete
with qGP.

Due to the computational cost of evaluations in ERL, only the 3x3 matrix
problems are considered. ERL on average finds the highest fitness solution at
a distinctly earlier rate than without obstacles. However, the strength of
this approach is found in the generalization ability of the evolved agents,
and Figure 4 simply plots the average best fitness achieved for a single
problem instance. The generalization ability of ERL agents is demonstrated in
Figure 5, which plots the number of unique belt problems solved by the single
best ERL agent at each generation (each line represents 1 of 10 independent
runs). The problems without obstacles are clearly easier, and the single best
agent generalizes to 10 unique problem configurations. When obstacles are
introduced, the best agent can only solve 4 unique problems.

(a) 3x3 without obstacles

(b) 3x3 with obstacles

Figure 5. Number of unique belt problems solved by the best ERL generalist
agent at each generation. A problem is considered solved if the agent manages
to connect the input and output with transport belts within 20 timesteps.

##  7\. Future Work

This paper represents an initial exploration into automating and finding high
quality solutions to in-game problems through an external interface. Our
future work will include a deeper analysis into the logistic belt problem,
expanding into the use of other logistic items including splitters, varying
speeds of logistics and updating our fitness functions to include
considerations for weighted constraints such as cost of materials, size of
solution, speed and throughput of transportation of materials.

Other problems in Factorio were considered at the beginning of our research.
One such problem is the electricity pole placement problem, which is
particularly interesting due to the relationship this problem has with
placement of machinery (as each item requires an electrical supply). This may
be an interesting problem to consider co-evolution techniques. Another example
of potential future work includes decision-making problems which consider
resource management, power management and pollution generation, and attempts
to control automated systems (such as the problem described in this paper) to
minimize or maximize certain objectives.

Each of the algorithms in this paper is tested upon the logistic belt problem.
Additional extensions to these algorithms could enhance their ability to cope
with large search spaces, such as path-finding hybridization. An important
feature of qGP is the ability to produce multiple outputs with a single
operator. Additional complex model structures may provide improvements to
tackling the logistic belt problem.

##  8\. Conclusions

The video game Factorio constitutes an untapped suite of research problems,
and this paper provides an initial step into exploring them. In this study we
present an introduction to problems which can be found within Factorio and
define one of these problems: the logistic belt optimization problem. We
provide three algorithms to attempt to solve the problem and compare their
efficacy. A mathematical model of the logistic belt optimization problem is
presented. This work also provides an interface with Factorio to allow
external optimization algorithms to connect to Factorio, allowing future
researchers to conduct experiments without experience with the in-game
scripting language, Lua. Our results show a variety of algorithm types and the
benefits of exploratory algorithms such as Parallel Simulated Annealing, the
strengths of domain-focused genetic programming and the flexibility of general
purpose metaheuristics, specifically evolutionary reinforcement learning.

##  9\. Acknowledgments

This work was supported in part by the National Institute of Food and
Agriculture (AFRI Project No. 2019-67015-29323), the BEACON Center for the
Study of Evolution in Action and through computational resources and services
provided by the Institute for Cyber-Enabled Research at Michigan State
University.

## References

  * (1) Source RCON Protocol, https://developer.valvesoftware.com/wiki/Source_RCON_Protocol
  * (2) Arafat, M.Y., Moh, S.: Routing protocols for unmanned aerial vehicle networks: A survey. IEEE Access 7, 99694–99720 (2019)
  * (3) Arzate Cruz, C., Ramirez Uresti, J.A.: Player-centered game AI from a flow perspective: Towards a better understanding of past trends and future directions. Entertainment Computing 20, 11–24 (2017)
  * (4) Asada, M., Stone, P., Veloso, M., Lee, D., Nardi, D.: Robocup: A treasure trove of rich diversity for research issues and interdisciplinary connections [tc spotlight]. IEEE Robotics Automation Magazine 26(3), 99–102 (2019)
  * (5) Blickle, T., Thiele, L.: A mathematical analysis of tournament selection. In: Eshelman, L. (ed.) Intl. Conference on Genetic Algorithms (ICGA-1995). pp. 9–15 (1995)
  * (6) Bloemhof-Ruwaard, J.M., Van Beek, P., Hordijk, L., Van Wassenhove, L.N.: Interactions between operational research and environmental management. European journal of operational research 85(2), 229–243 (1995)
  * (7) Boysen, N., Briskorn, D., Fedtke, S., Schmickerath, M.: Automated sortation conveyors: A survey from an operational research perspective. European Journal of Operational Research 276(3), 796–815 (2019)
  * (8) Cai, J.h., Meng, W.j., Chen, Z.m.: The optimisation research of screw conveyor. International Journal of Innovative Computing and Applications 3(3), 169–176 (2011)
  * (9) Carr, M.J., Wheeler, C.A., Robinson, P.W., Chen, B.: Reducing the energy intensity of overland conveying using a novel rail-running conveyor system. International Journal of Mining, Reclamation and Environment pp. 1–16 (2020)
  * (10) Chrabaszcz, P., Loshchilov, I., Hutter, F.: Back to basics: Benchmarking canonical evolution strategies for playing atari. CoRR abs/1802.08842 (2018), http://arxiv.org/abs/1802.08842
  * (11) Dantzig, G.B., Ramser, J.H.: The truck dispatching problem. Management science 6(1), 80–91 (1959)
  * (12) Dickinson, B.C., Jenkins, O.C., Moseley, M., Bloom, D., Hartmann, D.: Roomba pac-man: Teaching autonomous robotics through embodied gaming. In: Robots and Robot Venues: Resources for AI Education, Papers from the 2007 AAAI Spring Symposium, Technical Report SS-07-09, Stanford, California, USA, March 26-28, 2007. pp. 35–39. AAAI (2007), http://www.aaai.org/Library/Symposia/Spring/2007/ss07-09-008.php
  * (13) Duhan, S., Zhang, C., Jing, W., Li, M.: Factory optimization using deep reinforcement learning AI. Purdue Undergraduate Research Conference (2019), https://docs.lib.purdue.edu/purc/2019/Posters/57
  * (14) Fakhfakh, M., Cooren, Y., Sallem, A., Loulou, M., Siarry, P.: Analog circuit design optimization through the particle swarm optimization technique. Analog integrated circuits and signal processing 63(1), 71–82 (2010)
  * (15) Gavrilov, S., Zheleznikov, D., Khvatov, V., Chochaev, R.: Clustering optimization based on simulated annealing algorithm for reconfigurable systems-on-chip. In: 2018 IEEE Conference of Russian Young Researchers in Electrical and Electronic Engineering (EIConRus). pp. 1492–1495. IEEE (2018)
  * (16) Gerard, B., O’Rourke, L.: Optimisation of overland conveyor performance. Aust Bulk Handl Rev (January/February), 26–36 (2009)
  * (17) Gielen, G.G., Walscharts, H.C., Sansen, W.M.: Analog circuit design optimization based on symbolic simulation and simulated annealing. IEEE Journal of Solid-State Circuits 25(3), 707–713 (1990)
  * (18) Grinčová, A., MArAsová, D.: Experimental research and mathematical modelling as an effective tool of assessing failure of conveyor belts. Eksploatacja i Niezawodność 16(2) (2014)
  * (19) Jahir, Y., Atiquzzaman, M., Refai, H., Paranjothi, A., LoPresti, P.G.: Routing protocols and architecture for disaster area network: A survey. Ad Hoc Networks 82, 1–14 (2019)
  * (20) Kelly, S., Banzhaf, W.: Temporal memory sharing in visual reinforcement learning. In: Genetic Programming Theory and Practice XVII, pp. 0–0. Springer (2020)
  * (21) Kelly, S., Heywood, M.I.: Emergent solutions to high-dimensional multitask reinforcement learning. Evolutionary Computation 26(3), 378–380 (2018)
  * (22) Koza, J.R.: Genetic programming: on the programming of computers by means of natural selection. MIT Press 1 (1992)
  * (23) Kunze, T.: Video games and the education system – how they can benefit from each other. In: Elmenreich, W., Schallegger, R., Schniz, F., Gabriel, S., Pölsterl, G., Ruge, W. (eds.) Savegame, Perspektiven der Game Studies, pp. 31–40. Springer (2019)
  * (24) Leblanc, D.: The Only Game in Town: Simulators and the Circuits of Capitalism. Ph.D. thesis, Concordia University, Montreal, Canada (2018)
  * (25) Majerech, O.: Solving algorithms for multi-agent path planning with dynamic obstacles. Master’s thesis, Faculty of Mathematics and Physics, Charles University, Prague, Czech Republic (2017)
  * (26) Metropolis, N., Rosenbluth, A.W., Rosenbluth, M.N., Teller, A.H., Teller, E.: Equation of state calculations by fast computing machines. The Journal of Chemical Physics 21(6), 1087–1092 (1953)
  * (27) Mizutani, W.K., Kon, F.: Toward a reference architecture for economy mechanics in digital games. In: Proc. of the Brazilian Symposium on Games and Digital Entertainment (SBGames 2019) (2019)
  * (28) Mnih, V., Kavukcuoglu, K., Silver, D., Rusu, A.A., Veness, J., Bellemare, M.G., Graves, A., Riedmiller, M., Fidjeland, A.K., Ostrovski, G., Petersen, S., Beattie, C., Sadik, A., Antonoglou, I., King, H., Kumaran, D., Wierstra, D., Legg, S., Hassabis, D.: Human-level control through deep reinforcement learning. Nature 518(7540), 529–533 (2015)
  * (29) Phillips, A., Smith, G., Cook, M., Short, T.: Feminism and procedural content generation: toward a collaborative politics of computational creativity. Digital Creativity 27(1), 82–97 (2016)
  * (30) Pillac, V., Gendreau, M., Guéret, C., Medaglia, A.L.: A review of dynamic vehicle routing problems. European Journal of Operational Research 225(1), 1–11 (2013)
  * (31) Poli, R., McPhee, N.F., Vanneschi, L.: Elitism reduces bloat in genetic programming. In: Keijzer, M., Ryan, C., et al. (eds.) Proceedings of the 10th annual conference on Genetic and Evolutionary Computation - GECCO-2008. ACM Press, New York, USA (2008)
  * (32) Saha, B.K., Misra, S., Pal, S.: Seer: Simulated annealing-based routing in opportunistic mobile networks. IEEE Transactions on Mobile Computing 16(10), 2876–2888 (2017)
  * (33) Shaikh, F.S., Wismüller, R.: Routing in multi-hop cellular device-to-device (d2d) networks: A survey. IEEE Communications Surveys & Tutorials 20(4), 2622–2657 (2018)
  * (34) Silver, D., Schrittwieser, J., Simonyan, K., Antonoglou, I., Huang, A., Guez, A., Hubert, T., Baker, L., Lai, M., Bolton, A., Chen, Y., Lillicrap, T., Hui, F., Sifre, L., van den Driessche, G., Graepel, T., Hassabis, D.: Mastering the game of go without human knowledge 550, 354–359 (10 2017)
  * (35) Sutton, R.S., Barto, A.G.: Reinforcement Learning: An Introduction. MIT Press, 2nd edn. (2018)
  * (36) Vincent, F.Y., Redi, A.P., Hidayat, Y.A., Wibowo, O.J.: A simulated annealing heuristic for the hybrid vehicle routing problem. Applied Soft Computing 53, 119–132 (2017)
  * (37) Whalen, I., Banzhaf, W., Al Mamun, H.A., Gondro, C.: Evolving SNP Panels for Genomic Prediction. In: Banzhaf, W., Chen, B., Deb, K., Holekamp, K., Lenski, R., et al. (eds.) Evolution in Action: Past, Present and Future, pp. 467–487. Springer (2020)
  * (38) Yannakakis, G.N., Togelius, J.: Artificial Intelligence and Games. Springer (2018), http://gameaibook.org

  *[NPC]: Non-Player Character
  *[PCG]: Procedural Content Generation
  *[OR]: Operational Research
  *[ASS]: Automated Sortation Conveyor
  *[PSO]: Particle Swarm Optimization
  *[D2D]: Device to Device
  *[SC]: Soft Constraint
  *[FOI]: Factorio Optimizer Interface
  *[RCON]: Remote Console
  *[SA]: Simulated Annealing
  *[PSA]: Parallel Simulated Annealing
  *[GP]: Genetic Programming
  *[qGP]: q Genetic Programming
  *[RL]: Reinforcement Learning
  *[ERL]: Evolutionary Reinforcement Learning
  *[GA]: Genetic Algorithm
  *[TPG]: Tangled Program Graphs


# Highly tunable NbTiN Josephson junctions fabricated with focused helium ion
beam

Aki Ruhtinas<EMAIL_ADDRESS>Ilari J. Maasilta<EMAIL_ADDRESS>Nanoscience
Center, Department of Physics, University of Jyväskylä, FI-40014 Jyväskylä,
Finland

###### Abstract

We demonstrate a ”direct writing” method for the fabrication of planar
Josephson junctions from high quality superconducting niobium titanium nitride
(NbTiN) thin films using focused He-ion beam irradiation. Compared to the
materials previously used in such processing, YBCO and MgB2, NbTiN has much
better mechanical and electrical properties, as well as good corrosion
resistance. We show that we can control the suppression of superconductivity
in NbTiN as a function of the helium ion beam fluence, and that this
controllable critical temperature suppression combined with the high spatial
resolution and position control of the He-ion beam in a helium ion microscope
enables us to successfully fabricate Josephson junctions with highly tunable
weak links. Because of the continuous nature of the disorder-induced metal-
insulator transition, this method allows the creation of barriers with wide
range of resistivities ranging from the metallic to the insulating state, with
the critical current and the junction resistance varying over two orders of
magnitude. Electrical transport measurements show that junctions follow
closely the ideal resistively and capacitively shunted junction behavior, have
high characteristic voltages ($0.2-1.4$ mV) and show Shapiro steps up to very
high orders. This suggests that these type of junctions are suitable for a
wide range of applications in superconducting electronics and quantum
information technology, with the bonus that a whole device can be fabricated
from just a single thin film, with the excellent electrical and microwave
characteristics offered by NbTiN.

##  I Introduction

Weak disorder does not have an effect on the superconducting state [1], but
for strong enough disorder, with the mean free path and the Fermi-wavelength
of the same order of magnitude, a superconductor can undergo a superconductor-
to-insulator transition (SIT). This quantum phase transition can be exploited
in the fabrication of weak links for Josephson junctions, as changing the
level of disorder can change the film properties drastically. However, to make
a good junction, the disordered region needs to be defined with nanometer
resolution. Although ion irradiation or ion implantation of weak links is not
new idea [2], high enough spatial resolution became available with the
introduction of helium ion microscope (HIM), which became commercially
available only in 2007. This opened up new possibilities for nanofabrication,
as with HIM one can achieve a beam diameter of $\sim 0.5$ nm. In the direct-
write fabrication of Josephson junctions, this high resolution helium ion beam
is used to locally irradiate a superconducting film to introduce disorder and
therefore control the SIT, enabling fabrication of highly tunable weak links
where both the width and the strength of the link can be tuned continuously.
However, near the SIT, the superconducting film properties are relatively
sensitive to changes in the helium ion fluence, and therefore a very precise
fluence control is needed for an accurate tuning of the weak link.

In 2015, direct writing of Josephson junctions was first demonstrated in YBCO
films [3]. After this, the direct writing method has been applied not only to
YBCO thin films to create Josephson junctions and SQUIDS [4, 5, 6], but also
to MgB2 [7] and to Bi2Sr2CaCu2O8+x [8] thin films. However, these materials
are often difficult to work with, and electrical and mechanical properties are
not good enough for some applications. Also for many superconducting devices,
a smaller superconducting gap is desired. It is known that a disorder driven
superconductor-insulator transition existsin NbN [9, 10, 11], TiN [12, 13] and
NbTiN [14, 15], and thus direct writing should work in all of the materials,
in principle. Consequently, there has been demonstration of direct writing of
nanowires in NbN [16], and we have some preliminary experimental evidence that
direct writing of Josephson junctions should work in nitrides of NbN, NbTiN
and TiN [17].

Here, we focus on using NbTiN as the superconducting material, to which we
direct write Josephson junctions with a helium ion microscope beam. The
observed electrical characteristics of the junctions closely follow
established theory. We use pulsed laser deposition (PLD) to deposit high
quality NbTiN with a relatively high $T_{c}>15$ K. NbTiN has good electrical
and mechanical characteristics and high corrosion resistance which makes it an
excellent choice for Josephson junction fabrication. Moreover, NbTiN is
already a state-of-the art material in many devices, and it has been used in
microwave resonators [18, 19], rapid single flux quantum (RSFQ) devices [20],
superconducting nanowire single-photon detectors (SNSPDs)[21, 22] and in THz-
band heterodyne sensing technologies [23, 24]. As NbTiN is one of the most
promising materials in a wide variety of applications, the ability to
fabricate whole devices using only NbTiN seems to be a promising avenue to
pursue.

Figure 1: (a) a helium ion micrograph and (b) an atomic force microscope (AFM)
image of a Josephson junction created by helium ion irradiation. The fluence
used here was $\sim 6\cdot 10^{18}$ ion/cm2. The width of the irradiated
region is quite wide, $\sim$14 nm, and superconducting wire is partially
milled, as can be seen from both the AFM and HIM images. (c) Shows a SRIM
simulation of the film dislocations due to 30 kV helium ion beam impinging to
35 nm thick NbTiN2 target.

##  II Experimental details

###  II.1 Pulsed laser deposition

Thin films of NbTiN were deposited using the same reactive infrared pulsed
laser deposition (PLD) technique as demonstrated for TiN [25], TaN[26] and
NbN[27] before. We use a Q-switched Nd:YAG laser (Ekspla NL301 HT) operating
at its fundamental 1064 nm wavelength. The laser fluence at the target was
$\sim$6 J/cm2 and the pulse repetition frequency 10 Hz. As a target, we used
two 99.99% purity NbTi alloy targets (Matsurf Technologies Inc.) with two
different compositions, Nb0.39Ti0.61 and Nb0.5Ti0.5. The target was ablated in
ultra-high purity nitrogen atmosphere using $\sim 50$ mTorr pressure. The
substrate temperature was kept at $700\,^{\circ}$C, and the substrates were
glued to substrate heater using silver glue in order to maximize the thermal
contact. The growth rate of the film was $\sim$1 nm/min. As a substrate, we
used cubic (100)-oriented MgO (Crystec GmbH) for several reasons. First of
all, it has a good lattice match to NbTiN, allowing epitaxial growth. Second,
the SF6-based plasma etching step to define the wires stops to MgO. Third and
most importantly, MgO seems to handle large helium ion fluences without issues
with the formation of helium bubbles in the substrate, unlike Si, for example.
The superconducting critical temperature of our PLD deposited NbTiN films
depends on the thickness. for a 50 nm film thickness, we have obtained
$T_{c}\sim$15 K while the highest achieved $T_{c}$ is 15.7 K for a thickness
of 100 nm. For a 22 nm film, we have achieved a $T_{c}$ of 11.5 K. The
resistivity of the best films is $\sim$
$30\text{\,}\mathrm{\SIUnitSymbolMicro\SIUnitSymbolOhm}\text{\,}\mathrm{m}$,
and the superconducting transition width is typically quite narrow ($\sim$50
mK).

###  II.2 Electron beam lithography and plasma etching

Bonding pads and superconducting wires with widths ranging between
0.2-2$\text{\,}\mathrm{\SIUnitSymbolMicro m}$ were defined to the PLD
deposited nitride film with electron beam lithography and reactive ion
etching. Before the lithography, the chip was cleaned in isopropyl alcohol
(IPA) bath using an ultrasonic cleaner (Finnsonic) to get rid off the silver
glue residues and other contaminants. Then, a resist layer was spun on the
samples at 4000 rpm for 60 s, with the resulting resist layer thickness
approximately 400 nm. Here we used the negative e-beam resist AR-N 7520.17,
because it has excellent plasma etch resistance and because a negative resist
is more suitable for our device geometry. After spinning, an e-beam
lithography tool (Raith E-Line) was used to expose the resist to define the
pattern, with a
110-$150\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{C}\mathrm{/}\mathrm{c}\mathrm{m}^{2}$
dose depending on the pattern size. After the exposure, the resist was
developed in undiluted AZ 351B developer for 50s, and rinsed with IPA, and
reactive ion etching (RIE) was performed (Oxford Plasmalab 80+ RIE) to remove
the NbTiN film outside the wire and bonding pad regions, using 60 W RF power,
100,sccm SF6 flow, and 5 sccm O2 flow at a pressure of 70 mTorr. This results
in etching speed of roughly 15 nm/min. Finally, the remaining resist was
removed using O2 plasma cleaning in the same RIE tool.

Figure 2: Suppression of $T_{c}$ and the disorder driven SIT in NbTiN thin
films. Panel (a) shows the resistance of a superconducting NbTiN wire after
irradiation a region of the wire with a helium ion beam using different
fluences. The wire undergoes a transition from the superconducting to the
insulating phase as the helium ion fluence increases, with a suppression of
$T_{c}$ along the way. Panel (b) shows the $T_{c}$ of the irradiated region in
a semilogarithmic plot, showing a linear dependency of $T_{c}$ as a function
of the logarithm of the ion beam fluence. SIT region is shaded in blue. Panel
(c) shows a logarithmic dependence of the critical current density as a
function of the helium ion fluence.

###  II.3 Helium ion irradiation and simulation

Josephson junctions were defined to the superconducting wire using a helium
ion beam ”direct writing” technique. For this, we used $30\,$kV He+ beam in a
Zeiss Orion Nanofab helium ion microscope with the nominal spatial resolution
of 0.5 nm. We varied the beam parameters, but mostly we used a
$20\text{\,}\mathrm{\SIUnitSymbolMicro m}$ aperture and a spot size of 4,
resulting in a beam current of 0.5-1 pA. The junction was defined as a narrow
line irradiation, perpendicular to the length of the superconducting wire,
creating a weak link between the two electrodes (Figure 1). The irradiation
was performed with a 1 nm step size and typically a
$1000\text{\,}\mathrm{\SIUnitSymbolMicro s}$ dwell time for each step, with
the number of repeats adjusted to get the desired fluence. Depending on the
spot size and focusing, the width of the irradiated line varied between
$\sim$5 nm-500 nm. Larger irradiation widths were achieved by using large spot
size and by defocusing the beam significantly. Compared to the He-ion beam
written YBCO junctions [3], the nitride junctions in this study need roughly
two orders of magnitude more fluence to suppress superconductivity. Thus, in
this work a fluence in the range $10^{16}-10^{20}$ He+/cm2 was used.

After the irradiation, every junction was imaged with a quick HIM scan to
ensure that the junction was successfully defined and to check the resulting
linewidth. When correct HIM parameters are found, the linewidth is
reproducible between subsequent irradiations. Even though the ion beam current
of the instrument is not as stable as an electron beam in a SEM and fluctuates
over time, the time averaged current is stable enough for reproducible
junction properties.

In order to get an idea of the effects of the helium ion irradiation, we used
a simulation code for ion-matter interaction (SRIM) to simulate a $30\,$kV He+
beam impinging on a 35 nm NbTiN film on top of a MgO substrate (Figure 1).
These simulations show that most helium ions stop in the substrate (MgO) after
traveling roughly 120-200 nm, and that there is only small lateral straggle
when the ions travel through the NbTiN film. In Figure 1(c), we show the
dislocation density in NbTiN film caused by the helium ion beam, revealing
that dislocations are concentrated in very narrow region of a few nm, which is
good for an accurate junction definition.

However, the true value of the disorder is probably lower than what is shown
in Figure 1 (c). This is because at room temperature, substantial amount of
target damage is usually efficiently repaired by recrystallization.This effect
has been demonstrated with both experiments [28] and simulations [29], and it
is most prominent in metals, so a large deviation from SRIM simulations is
expected in our case. In addition, SRIM assumes a perfect crystal for each
impinging ion, and does not take into account sputtering and the accumulation
of dislocations. Because of the complex interplay between milling and creation
of disorder, we believe that only the dimensions of disordered region should
be determined from these simulations.

##  III Results

###  III.1 Suppression of superconductivity

In order to investigate the effects of helium ion irradiation on NbTiN thin
films first more generally, we have irradiated superconducting NbTiN wires
with HIM, using a varying fluence and size (linewidth) of the irradiated
region. For small linewidths ($\sim 5$ nm) we used a spot size 4, whereas for
large linewidths a spot size of 2 was used. Results are shown in Figure 2 (a),
revealing that helium ion beam irradiation induces strong enough local
disorder to suppresses superconductivity and ultimately push the irradiated
region through the superconductor-insulator transition. Superconducting
critical temperature suppression as a function of helium ion fluence $\Phi$ is
surprisingly linear in semilogarithmic scale [Fig. 2 (b)], indicating that
$T_{c}\propto\log(\Phi)$. This relation allows straightforward and accurate
$T_{c}$ control of the NbTiN film, which may be useful in different
applications. We should note that while it is difficult to determine the exact
width of a large area irradiation precisely (and thus the precise value of the
fluence), the width fluctuations between different irradiation rounds are
small, and a possible systematic error would only shift this linear
dependency.

While the irradiations with small widths ($\sim 5$ nm) have the property
$T_{c}\propto\log(\Phi)$ at lower fluences, the presence of Josephson
supercurrent makes it more difficult to determine $T_{c}$ of irradiated
region. However, we can use the critical supercurrent $I_{c}$ as a metric,
instead. With the small linewidth irradiations, there is sharp transition to
insulating phase beginning around $\sim 2\cdot 10^{18}$ ion/cm2 for 35 nm
thick films , and this transition shifts to a higher fluence ($\sim 1\cdot
10^{19}$ ion/cm2) for 97 nm films [Fig. 2 (b)]. This shift is expected, as
there are more atoms to displace when films are thicker. However, there is a
lot of energy left in He ions after traveling through the NbTiN, so that a
film thickness increase does not increase the transition fluence linearly. We
find that this SIT region (blue region in Figure 2 (b)) is the best for
Josephson junction fabrication, as in that case the critical current can be
tuned continuously from $\sim I_{c}$ of the bare wire down to zero. We
demonstrate this in Figure 2(c), where we plot the critical current density
$J_{c}$ of several junctions fabricated to $\sim 500$ nm wide and $35$ nm
thick wires, where the junctions are in this SIT region. We find that $J_{c}$
is exponentially proportional to the helium ion fluence, so that
$\log(I_{c})\propto\Phi$. Because the SIT region is relatively narrow in
$\Phi$, precise control of the helium ion beam is needed for reproducible
fabrication of junctions. This is especially true for the narrowest $<10$ nm
junctions, as small changes in the irradiation width change $\Phi$ and thus
$I_{c}$ drastically.

###  III.2 Josephson junctions

####  III.2.1 DC transport measurements

After we characterized the SIT as a function of fluence, we proceeded to
fabricate and electrically characterize weak links with varying strengths and
widths. These junctions were characterized via DC transport measurements down
to 1.5 K with a liquid helium evaporation cryostat. We measured the current-
voltage characteristics (IVC) at different temperatures using a battery
powered sweepbox and a voltage divider, to achieve a low noise level. Voltage
and current were measured using Ithaco preamplifiers (models 1201 and 1211),
and the differential conductance was simultaneously determined either with
lock-in technique or with numerical differentiation of the IVC. In Figure 3,
we show DC transport measurements of four representative NbTiN Josephson
junctions (JJ1-JJ4) fabricated with the direct writing method. In order to
extract more precisely the critical current ($I_{c}$) and normal state
resistance ($R_{n}$) values, the IVC of every junction was fitted at selected
temperatures with the RCSJ model [30]. From these fits, we extracted $I_{c}$,
$R_{n}$ and values for the junction capacitance $C$. The experimental
parameters of these junctions are presented in Table 1.

There are some uncertainties in the reported parameter values, as the junction
length and the milling depth are difficult to experimentally estimate
accurately. The length of the junction (width of the HIM irradiation region)
can be estimated from HIM and AFM images with reasonable accuracy and are
reported in Table 1, but the possible milling depth of the film was not
directly measured. As some milling of the film can evidently occur (an example
shown in Fig. 1(b)) , this can distort the reported junction area $A_{j}$ and
$J_{c}$ values in Table 1 (calculated assuming no milling), and these should
be regarded as upper and lower limits, respectively.

From the transport measurements we concluded that it is possible to tune the
critical current from $I_{c}$ of the superconducting leads continuously down
to zero, and thus $I_{c}$ can be controlled over several orders of magnitude.
In junctions JJ1-JJ4, $I_{c}$ varies from
$2.2\text{\,}\mathrm{\SIUnitSymbolMicro A}$(JJ1) to
$390\text{\,}\mathrm{\SIUnitSymbolMicro A}$(JJ4). The junction dimensions such
as the film thickness $t$, the superconducting wire width $W$ and the junction
length $L_{j}$ (=irradiation width) vary between junctions JJ1-JJ4, showing
flexibility of the method. JJ2 was fabricated from a 100 nm thick film,
demonstrating that direct writing method is feasible for such films, thicker
than the most commonly used $\sim$30 nm thickness, which to our knowledge has
not been reported before. In addition, we have demonstrated that fabrication
of longer junctions is also possible with the direct write method, as JJ3 is
30 nm long. For that case, the resistivity of the barrier is low, but the long
length limits the critical current.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

Figure 3: DC transport measurements of four Josephson junctions fabricated
using different helium ion fluences and different superconducting wire
geometries devices: JJ1 (a,b), JJ2 (c,d), JJ3 (e,f) and JJ4 (g,h). Left
column: current-voltage characteristics, right column: differential
conductance as a function of temperature and current. For all junctions,
comparison to different theories are shown in right column by the dashed
lines.

Calculated resistivities for the devices (barriers) ranged from $\sim 1\cdot
10^{-3}$$\text{\,}\mathrm{\SIUnitSymbolOhm}\text{\,}\mathrm{cm}$ (JJ4) to
$\sim 2\cdot 10^{-1}$$\text{\,}\mathrm{\SIUnitSymbolOhm}\text{\,}\mathrm{cm}$
(JJ1) (As milling of the junction was not taken into account, these are the
upper limit values). However if we assume an extreme case of the film
thickness milled down to 5 nm, we get only an order of magnitude smaller
resistivities, still high values for metals. Sheet resistances can be
calculated more accurately, and these vary between $\sim$0.3 k$\Omega$ (JJ4,
JJ3) to $\sim$50 k$\Omega$ (JJ1). In all of the junctions, $R_{n}$ seems to be
independent of temperature, indicating still metallic conductivity. In NbTiN,
a Bose metallic state has been experimentally detected when the sheet
resistance is above the superconducting resistance quantum $R_{q}$ [15], so
that $R_{\square}>h/4e^{2}=6.45$ k$\Omega$. When compared to the measured
normal metal $R_{\square}$ here, it seems that junctions JJ1 and JJ2 are in
the vicinity or above this limit, suggesting a possible presence of the Bose
metal phase. Sample JJ1 has high enough $R_{\square}$ that based on Ref. [14]
these films should already be in the insulating side of the SIT, indicating
that JJ1 could possibly be an overdamped SIS junction.

From the measured d$V$/d$I$(T,I) maps (Figure 3, right column) we can extract
the temperature dependent critical current $I_{c}(T)$ that can be used to
determine the barrier properties. When a barrier is strong, $I_{c}$ has weak
temperature dependence at low temperatures, and this gets weaker as the
barrier gets stronger. This is expected, as $I_{c}(T)$ scales exponentially in
case of transparent diffusive metallic barriers, and the saturation at low
temperatures is signature of either an insulating barrier, or a very strong
diffusive barrier [2]. We do not consider the clean limit here, as the barrier
is highly disordered due to the irradiation.

In order to get a more quantitative analysis of $I_{c}(T)$, we performed
numerical modeling, assuming diffusive transport through the weak link. The
relevant length scales for transport through the weak link are the mean free
path $\ell$ in the normal metal, and the normal metal coherence length
$\xi_{N}$, given by [2]

$\xi_{N}(T)=\sqrt{\frac{\hbar D}{2\pi k_{B}T}}$ (1)

where $D$ is the diffusivity. For junctions JJ1 and JJ2 we use the Kulik and
Omelyanchuk theory (KO-I) for short SNS junctions in the diffusive limit [31],
which applies when $\ell\ll L_{j}\ll\xi_{N}$. Within these assumptions, the
supercurrent is given by

$I(T,\varphi)=\frac{2\pi
k_{B}T}{eR_{n}}\sum_{\omega_{n}}\frac{2\Delta\cos(\varphi/2)}{\Omega_{n}}\arctan\frac{\Delta\sin(\varphi/2)}{\Omega_{n}},$
(2)

where $\Omega_{n}^{2}=(\hbar\omega_{n})^{2}+|\Delta|^{2}\cos^{2}(\varphi/2)$,
$\varphi$ is the phase difference, and Matsubara frequencies $\omega_{n}$ are
defined such that $\hbar\omega_{n}=\pi k_{B}T(2n+1)$. The superconducting gap
$\Delta$ was calculated with the BCS temperature dependence [30] and using
$\Delta(T=0)=2k_{B}T_{c}$ that is shown to be more suitable for nitrides.

To determine $I_{c}(T)$, we solved Eq. 2 numerically using $\varphi=1.25\pi/2$
and summed the Matsubara frequencies up to $n=10^{5}$. The results of these
calculations are shown in Figure 3 (right column) as dashed lines. For both
JJ1 and JJ2, the fits to KO-I theory are very good, suggesting that these
junctions are in the short diffusive limit, as one would expect from the
junction dimensions and the weak link strength. For JJ1, we also compare the
data to the simple Ambegaokar-Baratoff (AB) theory [32]

$I_{c}(T)=\frac{\pi\Delta(T)}{2eR_{n}}\tanh(\frac{\Delta(T)}{2k_{B}T})$ (3)

valid for good SIS tunnel junctions, but we observe that the fit (black dashed
line in Fig.3 (b)) is not as good as with the KO-I theory. From these
measurements we can conclude that most likely the weak link is behaving more
like a diffusive normal metal barrier than an insulating tunnel barrier. In
addition, the assumption $\ell\ll L_{j}\ll\xi_{N}$ seems to hold very well, as
the KO-I theory fits almost perfectly to the measured $I_{c}(T)$s. For JJ1,
the purely theoretical $I_{c}R_{n}$ value was significantly higher than the
measured one, and we therefore used the prefactor $A=0.14$ for the KO-I theory
and $A=0.18$ for the AB theory to get good fits.

For JJ2, the fitted current prefactor was $A=0.35$ for the KO-I theory, in
other words, the characteristic voltage was closer to the theoretical maximum.
While the KO-I theory fits the junction JJ2 data well, it does not give as
perfect a fit as for JJ1. This is most likely because the weak link in JJ2 is
not as perfectly in the short junction limit. This is even more true in the
case for junctions JJ3 and JJ4, where it is evident that now the approximation
for a short junction is not fulfilled. With junctions JJ3 and JJ4, the length
scales $L_{j},\xi_{N}$ are closer to each other, and the Thouless energy
$E_{th}$ defined as $E_{th}=\hbar D/L^{2}$ is closer to the superconducting
gap $\Delta$. In order to get theoretical predictions for $I_{c}$ in junctions
JJ3 and JJ4, a numerical simulation was performed using usadel1-code [33] that
solves the quasiclassical Usadel equations in the diffusive limit. Simulation
fits using this code are shown as dashed lines in Figures 3 (f) and (h). This
code was also used for junctions JJ1 and JJ2, with similar results to the KO-I
theory, as they are in the limit $\Delta/E_{th}\ll 1$.

The best fit to JJ3 $I_{c}(T)$ curve was achieved for $\Delta/E_{th}=5$, so
that $E_{th}\approx 250\,\mu$eV and $D\approx 3.4$ cm${}^{2}/$s. However, the
fit was not perfect. It was possible to reproduce either the behavior at low
temperatures or at high temperatures separately, but for the overall best fit
a compromise was needed, as shown in Fig. 3 (f). This is probably because with
$D$ determined from fit and using Eq. 1, we get $\xi_{N}(8\,\text{K})\approx$7
nm and $\xi_{N}(1.5\,\text{K})\approx$17 nm. These length scales begin to
approach the junction length (30 nm) at the low temperature end, causing most
likely a crossover between the short junction limit to the long junction limit
as a function of $T$. The correction factor for $I_{c}$ in the fit was
$A=0.21$, so the maximum characteristic voltage would be five times larger.
Interestingly, also for all other junctions fabricated with $W_{irr}\approx
30\,$nm and $I_{c}\sim 10-50\mu$A, we always have a similar $Ic(T)$ behavior
to JJ3, and we conclude that the $I_{c}(T)$ properties for JJ3 are common for
these type of junctions, although not very well described with the theory.

For the junction JJ4, we found the best fit for $\Delta/E_{th}=6$, so that we
have $E_{th}\approx 240\,\mu$eV, similar to the junction JJ3. As JJ4 junction
was shorter, the diffusion constant is only $D\approx 0.7$ cm${}^{2}/$s and
thus $\xi_{N}\sim 1-2$ nm at the temperatures of the experiment. This put
junction JJ4 well into the long junction limit. However, the calculated
current scaling factor was $A=1.8$, so that the measured current was almost a
factor of two higher than what was expected from the simulations. This
contribution is most likely due to excess current.

The characteristic voltage $I_{c}R_{n}$ values were also determined for each
junction. All junctions shown in Figure3 feature relatively high $I_{c}R_{n}$
values, ranging between $0.2-1.4$ mV. From the RCSJ fits can be seen that
excess current contribution in these junctions is mostly small, and for
example for junction JJ3, the excess current contribution is essentially zero.
This can been seen both from the RCSJ fits and the magnetotransport
measurements (Figure 5), where the critical current is completely suppressed
at the nodes of the Fraunhofer-pattern. The highest $I_{c}R_{n}$ junction is
JJ2, which has a characteristic voltage of $\sim$ 1.4 mV. One explanation for
the high value is that the sheet resistance of JJ2 is really close to $R_{q}$,
and there is evidence that $I_{c}R_{n}$ for weak links closer to the phase
transition are higher. Because of some nonlinearity in the IVC, it is
difficult to estimate the excess current precisely, but based on the RCSJ
fits, its contribution should be below 10%. Hence the high $I_{c}R_{n}$ is not
likely caused by an excess current contribution. The achieved $I_{c}R_{n}$
values are already suitable for example for RSFQ applications, and with some
tuning it should be possible to push these characteristic voltage values even
higher.

Because self-heating effects quite closely resemble capacitive effects [2,
34], the determination of capacitance from the RCSJ fits is known to be
inaccurate. We believe that the hysteresis seen in junction JJ4 is most likely
of thermal origin instead of capacitive. Nevertheless, thermal effects of
other junctions seem to be small, and when we take into account both the
geometric and the electronic capacitance $C_{e}=3\pi\hbar/(32\Delta
R_{n})$[35, 36], the measured values are roughly of the same order of
magnitude as the values determined from RCSJ fits. Although the above
expression for the electronic capacitance is derived for tunnel junctions,
this theory could be applicable here as well, as we have barriers in the
vicinity of the SIT.

Table 1: Junction parameters Name |  $T_{c}$(K) |  $t$ (nm) | w (nm) |  $\Phi$ (He+/cm2) | Lj(nm) | Aj (nm2) |  $I_{c}$ ($\mu$A) |  $R_{n}$ ($\Omega$) |  $I_{c}R_{n}$ (mV) |  $J_{c}$ (kA/cm2) | R${}_{\square}(\Omega$) | C(F)
---|---|---|---|---|---|---|---|---|---|---|---|---
JJ1 | 13.00 | 35 | 860 | 9.8$\cdot 10^{18}$ | 3 | 30100 | 2.2 | 173 | 0.381 | 7.31 | 29.8 | 1$\cdot 10^{-17}$
JJ2 | 15.70 | 100 | 300 | 1.3$\cdot 10^{19}$ | 12 | 30000 | 9.0 | 160 | 1.44 | 30.0 | 4.00 | 2$\cdot 10^{-15}$
JJ3 | 13.90 | 48 | 2030 | 1.1$\cdot 10^{19}$ | 30 | 97440 | 44 | 4.8 | 0.211 | 45.2 | 0.32 | 4$\cdot 10^{-13}$
JJ4 | 12.84 | 35 | 1500 | 2.9$\cdot 10^{18}$ | 14 | 52500 | 390 | 2.5 | 0.98 | 743 | 0.27 | 5$\cdot 10^{-14}$

Figure 4: RF characterization of junction JJ3, showing both 6 GHz (a,b) and 11
GHz (d,e) measurements and comparison to simulations (c,f). Panels (a) and (d)
show some selected IV curves with and without RF power, and theoretical fits
using the RSJ model. Panels (b) and (e) show the so called Shapiro map, i.e.
the measured $dI/dV$ as a function of both RF power and DC bias current. As a
comparison, full Shapiro map simulations are shown in (c) and (f). The best
fit parameters were found to be $R_{N}=4.8\,\Omega$, $Ic=17.5\,\mu$A.

Figure 5: Panel (a) shows the differential conductance of junction JJ3,
measured as a function of magnetic field at 4.2 K with a fit to Eq. 4(white
dotted line). There is clearly a periodic modulation of the critical current
as a function of the magnetic field, and at the nodes the critical current
vanishes completely. Panel (b) shows the same for junction JJ2, and a fit to
Eq.6. In this junction there is no such clear modulation as a function of
magnetic field.

####  III.2.2 Microwave measurements

In order to characterize the junctions further, we performed a
characterization of JJ3 with microwave irradiation. We used 2.1-11 GHz RF
radiation produced with Anritsu 68367C signal generator. The minimum signal
output was -10 dBm, but to get a lower power output for the measurement we
used a 20 dBm attenuator in the RF line. The RF signal was carried into the
cryostat via low loss coaxial cable, and the RF radiation was coupled to the
sample with simple loop antenna. The microwave measurements were carried out
by scanning the RF output power with a fixed frequency, and measuring $dV/dI$
directly using a four wire lock-in technique at each power level.

For a theoretical comparison, we performed RCSJ simulations for the IV curves,
and found that at 4.2 K junction behavior follows reasonably well the simpler
RSJ model, which was used to simulate the RF modulated response of the
junction vs. current and RF power (Shapiro map).

Figs. 4 (a) and (d) first shows examples of the IVCs of the junction JJ3 under
either 6 GHz or 11 GHz RF illumination, respectively, demonstrating clearly
Shapiro steps and agreement with the RSJ model simulations. The corresponding
full Shapiro maps are shown in Figs. 4 (b) and (e), showing nearly ideal
behavior. Comparison to the RSJ model simulated Shapiro maps (Fig. 4 (c) and
(f)) show very little deviation between the experiment and the simulation,
indicating a high quality junction. The Shapiro steps are visible up to high
order of $>20$ for the 6 GHz illumination.

####  III.2.3 Magnetic field measurements

Magnetotransport measurements were performed at 4.2 K using liquid helium
dipstick and a 4 T magnet. This measurement was performed by scanning the coil
current and measuring $dV/dI$ directly using a four wire lock-in technique.
For a rectangular junction geometry, the critical current is expected to
follow

$I_{c}(B_{z})=I_{c0}\left|\frac{\sin(\pi B_{z}A/\phi_{0})}{\pi
B_{z}A/\phi_{0}}\right|,$ (4)

where the junction area is $A=Lw$ and the magnetic field period for which the
phase changes by $2\pi$ is $\Delta B=\phi_{0}/A$. For the length of the
junction, the Josephson penetration depth $\lambda_{j}$ needs to be taken into
account, so that $L=L_{j}+2\lambda_{j}$ [30]. In Figure 5, we present the
results of these magnetotransport measurement for two junctions, JJ2 and JJ3.
Junction JJ3, made to a $2\text{\,}\mathrm{\SIUnitSymbolMicro m}$ wide
superconducting wire has clearly a sinusoidal modulation with a vanishing
critical current at the nodes. This proves that a Josephson current flows
through the junction, although the pattern slightly differs from the expected
$\sin(x)/x$ behaviour. Nodes are roughly equally spaced, but peculiarly, their
height does not decrease monotonically as the magnetic field increases. It is
known that while Equation 4 holds well for uniform current distribution and
for wide ($w\gg\lambda_{j}$) and thick films, in planar junctions made from
very thin films ($\lambda_{j}\gg t$) this diffraction pattern is often altered
because of flux focusing effects[37]. For very narrow junctions, the usual
magnetic period is replaced by $\Delta B\approx 1.84\Phi_{0}/w^{2}$, as was
shown by Rosenthal et al. [37]. These flux focusing effects can cause also
aperiodic node spacing, as well as nonmonotonically decreasing node height
[38], as observed in JJ3. For JJ3, the observed period is $\Delta B\approx
5.8$ mT, and if we calculate the Josephson penetration depth using this and
Eq. 4, we get 75$\pm$ 10 nm. This differs considerably from the usual
$\sim$200 nm reported in literature for NbTiN [39, 40]. The BCS prediction for
the London penetration depth is given by

$\lambda\approx
105\,\text{nm}\sqrt{\frac{\rho(\mu\Omega\text{cm})}{T_{c}(K)}},$ (5)

and using this we get $\lambda\approx 150$ nm for our best films, still a
factor of two larger. However, flux focusing effects increase the effective
length, and thus are not able to account for the discrepancy. Consequently,
for junction JJ3 the Rosenthal correction for flux focusing gives a way too
short node spacing ($\sim 1$ mT).

Magnetotransport measurements of junction JJ2, made from a 300 nm wide wire,
did not show the expected Franhofer diffraction pattern. The critical current
as a function of the applied magnetic field did not suppress completely, and
the shape of the $I_{c}(B)$ shows additional constant residual $I_{c}$ on top
of diffraction pattern (Fig.4 (b)). We surmise that this persisting $I_{c}$
may be caused by random structural inhomogenities to which Josephson vortices
are pinned [41, 30]. Yanson [41] showed that the shape of $I_{c}(H)$ is then
altered to

$I_{c}(H)=I_{c0}\sqrt{(1-\gamma^{2})\left(\sin(X)/X\right)^{2}+\gamma^{2}},$
(6)

where $X=\pi H/H_{0}$. Structural inhomogenities are introduced to this
equation by a random additional current $I_{1}(x)$ that is due to
inhomogenities in the barrier, and by a factor $N$ that represents the amount
of structural inhomogeneity so that
$\gamma^{2}=(\overline{I_{1}^{2}}/I_{0}^{2})\cdot(1/\pi N)$. When we fitted
Eq. 6 to the measured $I_{c}(B)$ of JJ2, we got a reasonable fit with
$\gamma=0.3$. From Figure 5 (b), we can determine the position of the first
minimum, and from this we get $\Delta B\approx 90$ mT. When we use Eq. 4 to
calculate the Josephson penetration depth for JJ2 we get 40$\pm$ 20 nm. Hence,
the determined $\lambda_{j}$ for both junctions are reasonably close to each
other, but differ considerably from the literature values as well as the BCS
prediction. For JJ2, the Rosenthal correction gives $\Delta B\approx 60$ mT,
which is not too far from the measured value and thus might explain the
observed node spacing in this case.

##  IV Conclusions

In this study, we have presented for the first time precise local control of
superconductivity of NbTiN films using He ion irradiation, and shown that this
can be utilized for the fabrication of highly tunable weak links for Josephson
junctions. We also extensively characterized the fabricated Josephson
junctions using both dc and microwave transport measurements, as well as
magnetotransport measurements, and show that the fabricated junctions are of
high quality. Additionally, we showed that we are able to tune the barrier
strength and fabricate really strong metallic barriers, pushing the limit
between metallic and insulating state. The junctions with strongest barrier
have excellent thermal stability (weak dependence on temperature) already at
1.5 K, which is quite promising for example for higher temperature quantum
information and superconducting electronics applications.

##  V Acknowledgements

###### Acknowledgements.

Authors wish to thank Geng Zhuoran and Pauli Virtanen for discussions and
providing parts of the codes, as well as Antti Kanniainen for the help with
the RF measurements. We also thank Tatu Korkiamäki for the help with some of
the cryogenic measurements. This research was supported by the Academy of
Finland project number 341823, the Vilho, Yrjö and Kalle Väisälä Foundation of
the Finnish Academy of Science and Letters and Jenny and Antti Wihuri
Foundation.

## References

  * Anderson [1959] P. Anderson, Theory of dirty superconductors, Journal of Physics and Chemistry of Solids 11, 26 (1959).
  * Likharev [1979] K. K. Likharev, Superconducting weak links, Rev. Mod. Phys. 51, 101 (1979).
  * Cybart  _et al._ [2015] S. Cybart, E. Cho, T. Wong, B. Wehlin, M. Ma, C. Huynh, and R. Dynes, Nano josephson superconducting tunnel junctions in YBa2Cu3O7–σ directly patterned with a focused helium ion beam, Nature nanotechnology 10, 598–602 (2015).
  * Müller  _et al._ [2019] B. Müller, M. Karrer, F. Limberger, M. Becker, B. Schröppel, C. Burkhardt, R. Kleiner, E. Goldobin, and D. Koelle, Josephson junctions and squids created by focused helium-ion-beam irradiation of ${\mathrm{y}\mathrm{ba}}_{2}{\mathrm{cu}}_{3}{\mathrm{o}}_{7}$, Phys. Rev. Applied 11, 044082 (2019).
  * Li  _et al._ [2020] H. Li, H. Cai, E. Y. Cho, S. J. McCoy, Y.-T. Wang, J. C. LeFebvre, Y. W. Zhou, and S. A. Cybart, High-transition-temperature nanoscale superconducting quantum interference devices directly written with a focused helium ion beam, Applied Physics Letters 116, 070601 (2020).
  * Chen  _et al._ [2022] Z. Chen, Y. Li, R. Zhu, J. Xu, T. Xu, D. Yin, X. Cai, Y. Wang, J. Lu, Y. Zhang, _et al._ , High temperature superconducting yba2cu3o7-$\delta$ josephson junction fabricated with a focused helium ion beam, Chinese Physics Letters (2022).
  * Kasaei  _et al._ [2018] L. Kasaei, T. Melbourne, V. Manichev, L. C. Feldman, T. Gustafsson, K. Chen, X. X. Xi, and B. A. Davidson, MgB2 josephson junctions produced by focused helium ion beam irradiation, AIP Advances 8, 075020 (2018).
  * Wang  _et al._ [2021] Y.-T. Wang, J. C. LeFebvre, E. Y. Cho, S. J. McCoy, H. Li, G. Gu, K. Kadowaki, and S. A. Cybart, Fabrication of Bi2Sr2CaCu2o8+x ab-plane josephson junctions by a focused helium ion beam, IEEE Transactions on Applied Superconductivity 31, 1 (2021).
  * Makise  _et al._ [2015] K. Makise, T. Odou, S. Ezaki, T. Asano, and B. Shinozaki, Superconductor–insulator transition in two-dimensional NbN/MgO and NbN/AlN/MgO films, Materials Research Express 2, 106001 (2015).
  * Ezaki  _et al._ [2012] S. Ezaki, K. Makise, B. Shinozaki, T. Odo, T. Asano, H. Terai, T. Yamashita, S. Miki, and Z. Wang, Localization and interaction effects in ultrathin epitaxial NbN superconducting films, Journal of Physics: Condensed Matter 24, 475702 (2012).
  * Yong  _et al._ [2013] J. Yong, T. R. Lemberger, L. Benfatto, K. Ilin, and M. Siegel, Robustness of the Berezinskii-Kosterlitz-Thouless transition in ultrathin NbN films near the superconductor-insulator transition, Phys. Rev. B 87, 184505 (2013).
  * Hadacek  _et al._ [2004] N. Hadacek, M. Sanquer, and J.-C. Villégier, Double reentrant superconductor-insulator transition in thin TiN films, Phys. Rev. B 69, 024505 (2004).
  * Baturina  _et al._ [2007] T. I. Baturina, C. Strunk, M. R. Baklanov, and A. Satta, Quantum Metallicity on the High-Field Side of the Superconductor-Insulator Transition, Phys. Rev. Lett. 98, 127003 (2007).
  * Burdastyh  _et al._ [2017] M. V. Burdastyh, S. V. Postolova, T. I. Baturina, T. Proslier, V. M. Vinokur, and A. Y. Mironova, Superconductor–Insulator Transition in NbTiN Films, JETP LETTERS 106, 749–753 (2017).
  * Burdastyh  _et al._ [2020] M. Burdastyh, S. Postolova, T. Proslier, S. Ustavshikov, A. Antonov, V. Vinokur, and A. Mironov, Superconducting phase transitions in disordered NbTiN films, Scientific Reports 10 (2020).
  * Martinez  _et al._ [2020] G. D. Martinez, D. Buckley, I. Charaev, A. Dane, D. E. Dow, and K. K. Berggren, Superconducting nanowire fabrication on niobium nitride using helium ion irradiation (2020).
  * Ruhtinas [2020] A. Ruhtinas,  _Direct writing of Josephson junctions in superconducting nitrides with focused helium ion beam_ , Master’s thesis, University of Jyväskylä (2020).
  * Bretz-Sullivan  _et al._ [2022] T. M. Bretz-Sullivan, R. M. Lewis, A. L. Lima-Sharma, D. Lidsky, C. M. Smyth, C. T. Harris, M. Venuti, S. Eley, and T.-M. Lu, High kinetic inductance nbtin superconducting transmission line resonators in the very thin film limit, Applied Physics Letters 121, 052602 (2022).
  * Müller  _et al._ [2022] M. Müller, T. Luschmann, A. Faltermeier, S. Weichselbaumer, L. Koch, G. B. P. Huber, H. W. Schumacher, N. Ubbelohde, D. Reifert, T. Scheller, F. Deppe, A. Marx, S. Filipp, M. Althammer, R. Gross, and H. Huebl, Magnetic field robust high quality factor nbtin superconducting microwave resonators, Materials for Quantum Technology 2, 015002 (2022).
  * Yu  _et al._ [2006] L. Yu, R. Gandikota, R. K. Singh, L. Gu, D. J. Smith, X. Meng, X. Zeng, T. V. Duzer, J. M. Rowell, and N. Newman, Internally shunted josephson junctions with barriers tuned near the metal–insulator transition for rsfq logic applications, Superconductor Science and Technology 19, 719 (2006).
  * Esmaeil Zadeh  _et al._ [2017] I. Esmaeil Zadeh, J. W. N. Los, R. B. M. Gourgues, V. Steinmetz, G. Bulgarini, S. M. Dobrovolskiy, V. Zwiller, and S. N. Dorenbos, Single-photon detectors combining high efficiency, high detection rates, and ultra-high timing resolution, APL Photonics 2, 111301 (2017).
  * Schuck  _et al._ [2013] C. Schuck, W. H. P. Pernice, and H. X. Tang, Nbtin superconducting nanowire detectors for visible and telecom wavelengths single photon counting on si3n4 photonic circuits, Applied Physics Letters 102, 051101 (2013).
  * Shiba  _et al._ [2012] S. Shiba, Y. Irimajiri, T. Yamakura, H. Maezawa, N. Sekine, I. Hosako, and S. Yamamoto, 3.1-thz heterodyne receiver using an nbtin hot-electron bolometer mixer and a quantum cascade laser, IEEE Transactions on Terahertz Science and Technology 2, 22 (2012).
  * Karpov  _et al._ [2007] A. Karpov, D. Miller, F. Rice, J. A. Stern, B. Bumble, H. G. LeDuc, and J. Zmuidzinas, Low noise 1 thz–1.4 thz mixers using nb/al-aln/nbtin sis junctions, IEEE Transactions on Applied Superconductivity 17, 343 (2007).
  * Torgovkin  _et al._ [2018] A. Torgovkin, S. Chaudhuri, A. Ruhtinas, M. Lahtinen, T. Sajavaara, and I. J. Maasilta, High quality superconducting titanium nitride thin film growth using infrared pulsed laser deposition, Superconductor Science and Technology 31, 055017 (2018).
  * Chaudhuri  _et al._ [2013] S. Chaudhuri, I. J. Maasilta, L. Chandernagor, M. Ging, and M. Lahtinen, Fabrication of superconducting tantalum nitride thin films using infrared pulsed laser deposition, Journal of Vacuum Science & Technology A 31, 061502 (2013), https://doi.org/10.1116/1.4812698 .
  * Chaudhuri  _et al._ [2011] S. Chaudhuri, M. R. Nevala, T. Hakkarainen, T. Niemi, and I. J. Maasilta, Infrared pulsed laser deposition of niobium nitride thin films, IEEE Transactions on Applied Superconductivity 21, 143 (2011).
  * Pelaz  _et al._ [2004] L. Pelaz, L. A. Marqués, and J. Barbolla, Ion-beam-induced amorphization and recrystallization in silicon, Journal of Applied Physics 96, 5947 (2004).
  * Nord  _et al._ [2002] J. Nord, K. Nordlund, and J. Keinonen, Amorphization mechanism and defect structures in ion-beam-amorphized si, ge, and gaas, Phys. Rev. B 65, 165329 (2002).
  * Tinkham [2004] M. Tinkham,  _Introduction to Superconductivity_ (Dover Publications, 2004).
  * Kulik and Omel’yanchuk [1975] I. O. Kulik and A. N. Omel’yanchuk, Contribution to the microscopic theory of the josephson effect in superconducting bridges, JETP Lett. (USSR) 21, 96 (1975).
  * Ambegaokar and Baratoff [1963] V. Ambegaokar and A. Baratoff, Tunneling between superconductors, Phys. Rev. Lett. 10, 486 (1963).
  * Virtanen and Heikkilä [2007] P. Virtanen and T. T. Heikkilä, Thermoelectric effects in superconducting proximity structures, Appl. Phys. A 89, 625 (2007), source code available at https://gitlab.jyu.fi/jyucmt/usadel1.
  * Courtois  _et al._ [2008] H. Courtois, M. Meschke, J. T. Peltonen, and J. P. Pekola, Origin of hysteresis in a proximity josephson junction, Phys. Rev. Lett. 101, 067002 (2008).
  * Eckern  _et al._ [1984] U. Eckern, G. Schön, and V. Ambegaokar, Quantum dynamics of a superconducting tunnel junction, Phys. Rev. B 30, 6419 (1984).
  * Larkin and Ovchinnikov [1983] A. I. Larkin and Y. N. Ovchinnikov, Decay of the supercurrent in tunnel junctions, Phys. Rev. B 28, 6281 (1983).
  * Rosenthal  _et al._ [1991] P. A. Rosenthal, M. R. Beasley, K. Char, M. S. Colclough, and G. Zaharchuk, Flux focusing effects in planar thin-film grain-boundary Josephson junctions, Applied Physics Letters 59, 3482 (1991).
  * Suominen  _et al._ [2017] H. J. Suominen, J. Danon, M. Kjaergaard, K. Flensberg, J. Shabani, C. J. Palmstrøm, F. Nichele, and C. M. Marcus, Anomalous fraunhofer interference in epitaxial superconductor-semiconductor josephson junctions, Phys. Rev. B 95, 035307 (2017).
  * Yu  _et al._ [2005] L. Yu, R. Singh, H. Liu, S. Wu, R. Hu, D. Durand, J. Bulman, J. Rowell, and N. Newman, Fabrication of niobium titanium nitride thin films with high superconducting transition temperatures and short penetration lengths, Applied Superconductivity, IEEE Transactions on 15, 44 (2005).
  * Yamamori  _et al._ [2010] H. Yamamori, H. Sasaki, and S. Kohjiro, Preparation of overdamped NbTiN Josephson junctions with bilayered Ti–TiN barriers, Journal of Applied Physics 108, 113904 (2010).
  * Rosenthal  _et al._ [1970] P. A. Rosenthal, M. R. Beasley, K. Char, M. S. Colclough, and G. Zaharchuk, Effect of fluctuations on the dependence of the Josephson current on magnetic field strength, Zh. Eksp. Teor. Fiz 58, 1497 (1970).