# HOBOTAN: Efficient Higher Order Binary Optimization Solver with Tensor
Networks and PyTorch

Shoya Yasuda<EMAIL_ADDRESS>Vignette & Clarity, Inc., 6-24-2
Honkomagome, Bunkyo-ku, Tokyo, 113-0021, Japan  Tokyo Institute of Technology,
School of Computing, J2 bldg. room1710, 4259 Nagatsuta-cho, Midori-ku,
Yokohama, 226-8503, Japan  Shunsuke Sotobayashi<EMAIL_ADDRESS>Independent Researcher  Yuichiro Minato<EMAIL_ADDRESS>blueqat Inc.,
2-24-12 Shibuya, Shibuya-ku, Tokyo 150-6139, Japan

###### Abstract

In this study, we introduce HOBOTAN, a new solver designed for Higher Order
Binary Optimization (HOBO). HOBOTAN supports both CPU and GPU, with the GPU
version developed based on PyTorch, offering a fast and scalable system. This
solver utilizes tensor networks to solve combinatorial optimization problems,
employing a HOBO tensor that maps the problem and performs tensor contractions
as needed. Additionally, by combining techniques such as batch processing for
tensor optimization and binary-based integer encoding, we significantly
enhance the efficiency of combinatorial optimization. In the future, the
utilization of increased GPU numbers is expected to harness greater
computational power, enabling efficient collaboration between multiple GPUs
for high scalability. Moreover, HOBOTAN is designed within the framework of
quantum computing, thus providing insights for future quantum computer
applications. This paper details the design, implementation, performance
evaluation, and scalability of HOBOTAN, demonstrating its effectiveness.

††preprint: APS/123-QED

##  I Introduction

Combinatorial optimization problems play a crucial role in various fields such
as logistics [1], finance, healthcare, and energy management. Solving these
problems requires advanced computational resources [2] and efficient
algorithms.

In this study, we developed a new solver called HOBOTAN, which supports Higher
Order Binary Optimization (HOBO). HOBOTAN is designed to solve high-order
combinatorial optimization problems using tensor networks [3] and supports
both CPU and GPU. The GPU version, based on PyTorch, is particularly designed
to be fast and scalable, showing high performance even for large optimization
problems.

The core feature of HOBOTAN is the preparation of a HOBO tensor that maps the
problem and performs tensor contractions as needed. This enables efficient
solving of complex combinatorial optimization problems. Furthermore, by
combining batch processing for tensor optimization and binary-based integer
encoding techniques, we significantly extend the advantages of combinatorial
optimization.

In the future, increasing the number of GPUs is expected to harness greater
computational power, allowing for high scalability by efficiently coordinating
multiple GPUs. This will enable HOBOTAN to handle even larger optimization
problems beyond the current capabilities of high-performance computing.

Additionally, HOBOTAN is designed within the framework of quantum computing,
providing insights for future quantum computer applications. When quantum
computers become practical, the optimization techniques using tensor networks
developed in HOBOTAN can potentially be applied directly. The knowledge gained
from the development of HOBOTAN is expected to contribute to research in the
field of quantum computing.

This paper first details the design and implementation of HOBOTAN, followed by
performance evaluation through some benchmark tests. Finally, we discuss the
new possibilities HOBOTAN brings to combinatorial optimization and outline
future research challenges and prospects.

##  II Theory

###  II.1 Overview of Tensor Networks

Tensor networks [4] are a mathematical framework designed to efficiently
handle complex operations on high-dimensional tensors. A tensor is a multi-
dimensional array, generalizing scalars (0-dimensional), vectors
(1-dimensional), and matrices (2-dimensional). Tensor networks reduce
computational complexity by contracting tensor elements. Contraction is an
operation that sums over specific dimensions of tensors, thus aggregating into
smaller tensors while maintaining the overall relationships within the
network.

###  II.2 Formulation of QUBO and HOBO in Combinatorial Optimization Problems

Combinatorial optimization problems involve finding a combination of variables
that minimizes (or maximizes) a given cost function [5]. Quadratic
Unconstrained Binary Optimization (QUBO) is formulated as follows:

$f(x)=\sum_{i}a_{i}x_{i}+\sum_{i<j}b_{ij}x_{i}x_{j}$

where $x_{i}$ are binary variables (0 or 1), and $a_{i}$ and $b_{ij}$ are
coefficients.

Higher Order Binary Optimization (HOBO) generalizes QUBO by including higher-
order terms. HOBO is formulated as follows:

$f(x)=\sum_{i}a_{i}x_{i}+\sum_{i<j}b_{ij}x_{i}x_{j}+\sum_{i<j<k}c_{ijk}x_{i}x_{j}x_{k}+\cdots$

where $c_{ijk}$ are third-order coefficients, and higher-order terms are
included.

###  II.3 Formulating HOBO as Multi-dimensional Tensors and Implementing
Contractions with Vector x

Based on the HOBO formulation, complex combinatorial optimization problems
with higher-order terms can be implemented using tensor networks. First,
prepare the tensor (HOBO tensor) representing the HOBO problem, where each
element of the tensor corresponds to a coefficient of the problem.

Next, define the binary variables $x$ as a vector and perform contractions
with the HOBO tensor. Contraction sums over the dimensions of the tensor and
the elements of the vector, calculating the cost function value.

Specifically, this can be formulated as:

$H(x)=\sum_{i,j,k,\ldots}H_{ijk\ldots}x_{i}x_{j}x_{k}\ldots$

where $H_{ijk\ldots}$ are the elements of the HOBO tensor, and
$x_{i},x_{j},x_{k},\ldots$ are binary variables. By contracting tensor $H$
with vector $x$, the cost function value is computed.

The implementation steps are as follows:

  1. 1.

Prepare the tensor $H$ representing the problem.

  2. 2.

Define the binary variable vector $x$.

  3. 3.

Contract tensor $H$ with vector $x$ to compute the cost function value.

###  II.4 Methods for Minimizing or Maximizing the Cost Function

After computing the cost function value, various optimization algorithms can
be used to minimize or maximize the cost function, thus solving the
combinatorial optimization problem. The following are some representative
methods:

####  II.4.1 Simulated Annealing

Simulated annealing is a probabilistic optimization algorithm inspired by the
annealing process in physics. It starts at a high temperature and gradually
cools down, reducing the search space. In the initial stages, a wide search
prevents getting stuck in local optima, eventually approaching the global
optimum. Simulated annealing can efficiently find optimal solutions from a
broad solution space for HOBO problems.

####  II.4.2 Gradient Descent

Gradient descent is an optimization method that uses the gradient (slope) of
the cost function to minimize its value. The gradient of the cost function is
computed, and variables are updated along the gradient direction to converge
to the optimum. For HOBO problems, the gradient is computed based on tensor
contraction results, and variables are updated iteratively. Variants like
mini-batch gradient descent can improve computational efficiency for large-
scale problems.

These methods enable efficient optimization of HOBO problems using tensor
networks. Thus, HOBOTAN serves as a powerful tool for solving complex
combinatorial optimization problems.

##  III HOBOTAN Implementation

###  III.1 Overview of Einsum Calculation

Einsum (Einstein summation convention) is a powerful notation used for tensor
operations, particularly in tensor network calculations. It allows for concise
and efficient representation of complex tensor contractions, which are crucial
for solving high-dimensional optimization problems. Einsum simplifies the
expression of these operations by summing over repeated indices in a product
of tensors.

For example, the matrix multiplication $C=AB$ can be expressed using Einsum as
follows:

$C_{ij}=\sum_{k}A_{ik}B_{kj}$

This efficiently performs tensor contractions, simplifying the description of
tensor calculations and improving computational efficiency.

###  III.2 Utilizing PyTorch for Combinatorial Optimization

Considering future scalability, especially leveraging GPU capabilities, we
utilize PyTorch, a machine learning framework, for combinatorial optimization
calculations. By using PyTorch as the backend, we benefit from parallel
computation and automatic differentiation, which are advantageous for solving
optimization problems. In PyTorch, tensors are created using torch.tensor, and
Einsum calculations are performed using torch.einsum.

###  III.3 Conversion to HOBO Tensor

To solve HOBO problems, the problem must be converted into a HOBO tensor.
Specifically, combinatorial optimization problems with high-order terms are
transformed into tensors of appropriate dimensions. In this conversion, it is
crucial to replicate the smallest subscript number to extend the dimensions.

For example, converting a 3-dimensional tensor term $c_{ijk}x_{i}x_{j}x_{k}$
into a 5-dimensional HOBO tensor involves replicating the smallest subscript
number to extend the dimensions as follows:

$c_{ijk}x_{i}x_{j}x_{k}\rightarrow c_{iiijk}x_{i}x_{i}x_{i}x_{j}x_{k}$

##  IV Special Implementations in HOBOTAN

###  IV.1 Compilation of Formulation

As mentioned in the previous chapter, if the order of terms does not reach the
maximum order of the HOBO tensor, the compiler automatically replicates the
smallest order to store the coefficients in the HOBO tensor. Specifically,
when converting a 3-dimensional term $c_{ijk}x_{i}x_{j}x_{k}$ into a
5-dimensional HOBO tensor, the smallest subscript number is replicated to
extend the dimensions. This results in the following transformation:

$c_{ijk}x_{i}x_{j}x_{k}\rightarrow c_{iiijk}x_{i}x_{i}x_{i}x_{j}x_{k}$

This transformation ensures that the HOBO tensor has consistent dimensions and
can appropriately represent all terms.

###  IV.2 Integer Encoding

Using HOBO allows for integer encoding from binary encoding, significantly
reducing the number of qubits required for the algorithm. This method defines
a new integer variable $y_{i}$ using the original binary variables $x_{i}$.
Specifically, this is expressed as:

$y_{i}=\sum_{k=0}^{n}2^{k}x_{ik}$

where $x_{ik}$ are binary variables and $y_{i}$ takes integer values. This
encoding enhances the expressive power of the problem and reduces the required
number of qubits.

###  IV.3 Tensor Batch System

For cost calculation in HOBO (including QUBO), we use tensor contraction
through einsum. Typically, tensor contraction calculates the cost as a scalar
by contracting with the corresponding vector $x$. In this implementation, to
fully utilize GPU performance and achieve parallelization, we prepare multiple
vectors $x$, expand them into a matrix, and contract them directly with the
HOBO tensor. For example, if the matrix $X$ contains multiple vectors $x$, the
cost can be calculated as follows:

$E=\sum_{i,j,k,\ldots}H_{ijk\ldots}X_{i}X_{j}X_{k}\ldots$

This method allows us to simultaneously compute the cost corresponding to
multiple vectors. This approach supports multiple sampling typically executed
by heuristic algorithms, ensuring parallelism and speeding up computation
through tensor network-based techniques.

###  IV.4 Tensor Dimension Reduction and Efficiency Through Tensor Networks

The core HOBO tensor of this algorithm is computationally intensive. By
optimizing this configuration, we can significantly reduce the computational
load. Specifically, unnecessary dimensions are reduced, and tensor network
optimization is performed to improve computational efficiency.

##  V Solving Specific HOBO Problems

Next, we demonstrate how to solve specific HOBO problems using the solver.
Three examples are provided here.

###  V.1 Problem 1: Maximizing the Number of Students Seated in a 5 by 5 Grid

In this problem, students cannot be seated in three consecutive seats either
in a row or in a column, which is a 3rd-order constraint.

####  V.1.1 Implementation

⬇

1import numpy as np

2from hobotan import *

3import matplotlib.pyplot as plt

4

5# Define quantum bits

6q = symbols_list([5, 5], ’q{}_{}’)

7

8# Objective: Preferably occupy all seats

9H1 = 0

10for i in range(5):

11 for j in range(5):

12 H1 += - q[i, j]

13

14# Constraint: No three consecutive seats in any row or column

15H2 = 0

16for i in range(5):

17 for j in range(5 - 3 + 1):

18 H2 += np.prod(q[i, j:j+3])

19for j in range(5):

20 for i in range(5 - 3 + 1):

21 H2 += np.prod(q[i:i+3, j])

22

23# Combine terms

24H = H1 + 10 * H2

25

26# Compile into HOBO tensor

27hobo, offset = Compile(H).get_hobo()

28print(f’offset\n{offset}’)

29

30# Select sampler

31solver = sampler.MIKASAmpler()

32

33# Perform sampling

34result = solver.run(hobo, shots=10000)

35

36# Display top 3 results

37for r in result[:3]:

38 print(f’Energy {r[1]}, Occurrence {r[2]}’)

39

40 # Convert result to array

41 arr, subs = Auto_array(r[0]).get_ndarray(’q{}_{}’)

42 print(arr)

43

44 # Visualize result

45 img, subs = Auto_array(r[0]).get_image(’q{}_{}’)

46 plt.figure(figsize=(2, 2))

47 plt.imshow(img)

48 plt.show()

####  V.1.2 Results

    offset
    0
    MODE: GPU
    DEVICE: cuda:0
    Energy -17.0, Occurrence 844
    [[1 1 0 1 1]
     [1 1 0 1 1]
     [0 0 1 0 0]
     [1 1 0 1 1]
     [1 1 0 1 1]]
    Energy -17.0, Occurrence 498
    [[1 1 0 1 1]
     [0 1 1 0 1]
     [1 0 1 1 0]
     [1 1 0 1 1]
     [0 1 1 0 1]]
    Energy -17.0, Occurrence 538
    [[1 1 0 1 1]
     [1 0 1 1 0]
     [0 1 1 0 1]
     [1 1 0 1 1]
     [1 0 1 1 0]]

The solver successfully identified configurations that maximize the number of
students seated while adhering to the constraint of not having three
consecutive seats occupied in any row or column. By introducing the tensor
batch system, the solver efficiently handles real-world problems, allowing it
to solve even 10,000 samples rapidly. This demonstrates the effectiveness and
utility of the HOBO solver in tackling combinatorial optimization problems
with higher-order constraints.

Figure 1: Configuration 1: Energy -17.0, Occurrence 844 Figure 2:
Configuration 2: Energy -17.0, Occurrence 498 Figure 3: Configuration 3:
Energy -17.0, Occurrence 538

###  V.2 Problem 2: Finding Pythagorean Triples

Next, we solve the problem of finding all sets of three integers that satisfy
the relationship $a^{2}+b^{2}=c^{2}$, known as Pythagorean triples.

####  V.2.1 Implementation

⬇

1import numpy as np

2from hobotan import *

3

4# Define quantum bits

5q = symbols_list([3, 4], ’q{}_{}’)

6

7# Integer Encoding

8x = 0

9y = 0

10z = 0

11for i in range(4):

12 x += (2**i) * q[0][i]

13 y += (2**i) * q[1][i]

14 z += (2**i) * q[2][i]

15

16# Objective : Pythagoras cost

17H1 = (x**2 + y**2 - z**2) **2

18

19# Constraint: For the qubits composing x, none of the qubits are set to 0.
The same applies to y and z

20H2 = 0

21H2 += np.prod((1 - q[0, :]))

22H2 += np.prod((1 - q[1, :]))

23H2 += np.prod((1 - q[2, :]))

24

25# Combine terms

26H = H1 + 10*H2

27

28# Compile into HOBO tensor

29hobo, offset = Compile(H).get_hobo()

30print(f’offset\n{offset}’)

31

32# Select sampler

33solver = sampler.MIKASAmpler()

34

35# Perform sampling

36result = solver.run(hobo, shots=10000)

37

38# Display top 5 results

39for r in result[:5]:

40 print(f’Energy {r[1]}, Occurrence {r[2]}’)

41

42 # Display the result for each integer

43 print("x =", Auto_array(r[0]).get_nbit_value(x))

44 print("y =", Auto_array(r[0]).get_nbit_value(y))

45 print("z =", Auto_array(r[0]).get_nbit_value(z))

####  V.2.2 Results

    offset
    30.0
    MODE: GPU
    DEVICE: cuda:0
    Energy -30.0, Occurrence 480
    x = 8.0
    y = 6.0
    z = 10.0
    Energy -30.0, Occurrence 842
    x = 4.0
    y = 3.0
    z = 5.0
    Energy -30.0, Occurrence 1430
    x = 12.0
    y = 9.0
    z = 15.0
    Energy -30.0, Occurrence 1541
    x = 12.0
    y = 5.0
    z = 13.0
    Energy -30.0, Occurrence 150
    x = 6.0
    y = 8.0
    z = 10.0

The program above defines a function to create custom binary variables and
constraints to ensure non-zero values. It then sets up the HOBO problem to
find Pythagorean triples, compiles it into a HOBO tensor, and uses the
MIKASAmpler to find solutions. The results show the first 8 solutions found,
along with their corresponding energies and occurrences, and the execution
time.

###  V.3 Problem 3: Solving the Traveling Salesman Problem using Integer
Encoding

####  V.3.1 Problem Formulation

In this problem, we solve the Traveling Salesman Problem (TSP) using integer
encoding to represent the path and distance [6]. The objective is to find the
shortest possible route that visits each city once and returns to the origin
city. While binary encoding typically uses one-hot constraints, integer
encoding has been shown to significantly reduce the number of required qubits,
allowing for more efficient solutions.

In the traditional TSP, the order of visiting cities is represented by one-hot
vectors. However, with integer encoding, the sequence of visiting cities is
represented as integers from 0 to 3. For this implementation, we set
constraints such that the product of the integers representing the sequence of
the 1st, 2nd, and 3rd cities equals 6, following the constraints from existing
literature.

####  V.3.2 Implementation

⬇

1import numpy as np

2from hobotan import *

3

4# Define quantum bits

5q = symbols_list([4, 2], ’q{}_{}’)

6xB = 2*q[0][0] + q[0][1]

7xC = 2*q[1][0] + q[1][1]

8xD = 2*q[2][0] + q[2][1]

9

10H1 = 0 # This time, we will not set a cost

11

12H2 = (xB*xC*xD-6)**2

13

14# Combine terms

15H = H1 + 10*H2

16

17# Compile into HOBO tensor

18hobo, offset = Compile(H).get_hobo()

19print(f’offset\n{offset}’)

20

21# Select sampler

22solver = sampler.MIKASAmpler()

23

24# Perform sampling

25result = solver.run(hobo, shots=10000)

26

27# Display top 3 results

28for r in result[:3]:

29 print(f’Energy {r[1]}, Occurrence {r[2]}’)

30

31 # Display the result for each number

32 print("xB =", Auto_array(r[0]).get_nbit_value(xB))

33 print("xC =", Auto_array(r[0]).get_nbit_value(xC))

34 print("xD =", Auto_array(r[0]).get_nbit_value(xD))

    offset
    360.0
    MODE: GPU
    DEVICE: cuda:0
    Energy -360.0, Occurrence 865
    xB = 1.0
    xC = 2.0
    xD = 3.0
    Energy -360.0, Occurrence 1886
    xB = 1.0
    xC = 3.0
    xD = 2.0
    Energy -360.0, Occurrence 780
    xB = 2.0
    xC = 1.0
    xD = 3.0

We were able to solve the TSP efficiently by using integer encoding for the
city visit sequences, significantly reducing the number of required qubits and
improving computational efficiency.

##  VI Optimizing HOBO Tensor

In this section, we discuss the optimization of the HOBO tensor and the use of
tensor networks for contraction and decomposition. The optimization of large
tensors is crucial because they can be computationally intensive. By carefully
designing the contraction path, we can achieve faster computations.

###  VI.1 Tensor Network Contraction

Tensor networks are graphical representations of tensor contractions. Each
node represents a tensor, and each edge represents a contraction between two
tensors. The goal is to contract all tensors to obtain a final scalar or
reduced tensor.

For example, consider the contraction of three tensors $A$, $B$, and $C$:

$D=\sum_{i,j,k}A_{ij}B_{jk}C_{ki}$

In this case, the order of contractions significantly affects the
computational complexity. We can optimize the contraction path to minimize the
total number of operations.

###  VI.2 Optimizing Contraction Paths

The choice of contraction path is crucial for efficient tensor computations.
Different network shapes and configurations require different optimal paths.
Tools like the opt_einsum library in Python help find the best contraction
path automatically.

For example, in the contraction of a high-dimensional tensor network, the
optimal path can be chosen based on the following criteria:

  * •

Minimize the intermediate tensor size: Reduce the size of tensors generated
during intermediate steps.

  * •

Reduce the total number of operations: Minimize the computational cost by
choosing the path with the least number of contractions.

By leveraging these optimization techniques, we can significantly improve the
performance of the HOBO solver, making it capable of handling larger and more
complex problems efficiently.

###  VI.3 Decomposing Large Tensors and Preparing for Path Optimization

By decomposing large tensors, we can increase the options for selecting tensor
contraction paths. Additionally, introducing approximations during
decomposition can further reduce computational weight. In this section, we
prepare for path optimization through decomposition without approximation. The
primary decomposition method we will use is Singular Value Decomposition
(SVD).

###  VI.4 Singular Value Decomposition (SVD)

Singular Value Decomposition (SVD) is a method for decomposing any matrix $A$
into the product of three simpler matrices:

$A=U\Sigma V^{*}$

Where:

  * •

$A$ is an $m\times n$ matrix.

  * •

$U$ is an $m\times m$ unitary matrix.

  * •

$\Sigma$ is an $m\times n$ diagonal matrix with singular values on the
diagonal.

  * •

$V^{*}$ is the conjugate transpose of an $n\times n$ unitary matrix.

This decomposition expresses the original matrix $A$ as a product of three
simpler matrices, simplifying tensor operations.

###  VI.5 Tensor Train Decomposition (TT Decomposition)

Next, we use SVD consecutively to decompose a large tensor into a network of
smaller tensors, known as Tensor Train Decomposition (TT Decomposition). TT
Decomposition expresses a high-dimensional tensor as a sequence of
3-dimensional tensors.

When decomposing a high-dimensional tensor $\mathcal{A}$ into a TT format, it
can be represented as:

$\mathcal{A}_{i_{1},i_{2},\ldots,i_{d}}=\sum_{r_{0},r_{1},\ldots,r_{d-1}}G_{1}[r_{0},i_{1},r_{1}]\cdots
G_{d}[r_{d-1},i_{d},r_{d}]$

Where:

  * •

$\mathcal{A}$ is a $d$-dimensional tensor.

  * •

$G_{k}$ is a 3-dimensional tensor such that $G_{k}\in\mathbb{R}^{r_{k-1}\times
n_{k}\times r_{k}}$.

  * •

$r_{0}$ and $r_{d}$ are virtual boundary conditions set to 1.

This TT Decomposition is performed sequentially as follows:

  1. 1.

Reshape the high-dimensional tensor $\mathcal{A}$ into a matrix $A_{(1)}$.

  2. 2.

Apply SVD to $A_{(1)}$, decomposing it into $U\Sigma V^{*}$.

  3. 3.

Incorporate $\Sigma V^{*}$ into the next tensor and reshape it into a matrix
$A_{(2)}$.

  4. 4.

Repeat this process until the $d$-dimensional tensor is fully decomposed.

Specifically, the process can be written as:

$A_{(1)}\approx U_{1}\Sigma_{1}V_{1}^{*}$ $A_{(2)}=\Sigma_{1}V_{1}^{*}$
$A_{(2)}\approx U_{2}\Sigma_{2}V_{2}^{*}$

By repeating this process, a large tensor is decomposed into a network of
smaller tensors. This approach increases the options for contraction paths and
lays the groundwork for optimization.

###  VI.6 Solving TSP Using TT Decomposition

In this subsection, we explore solving the Traveling Salesman Problem (TSP)
using Tensor Train (TT) decomposition. By applying TT decomposition, we can
significantly reduce the computational complexity of the problem.

####  VI.6.1 Before Decomposition

In the TSP problem, we used the same formulation as in the previous section’s
Problem 3: Solving the Traveling Salesman Problem using Integer Encoding. The
size of the HOBO tensor was $(6,6,6,6,6,6)$.

This is the result of the numpy einsum path when contracting this tensor
directly with the vector $x$ to obtain the answer:

      Complete contraction:  i,j,k,l,m,n,ijklmn->
             Naive scaling:  6
         Optimized scaling:  6
          Naive FLOP count:  3.266e+05
      Optimized FLOP count:  1.120e+05
       Theoretical speedup:  2.917
      Largest intermediate:  7.776e+03 elements

  * •

Complete contraction: All indices $i,j,k,l,m,n$ are contracted to obtain the
final scalar value.

  * •

Naive scaling: The scaling of the calculation using a naive approach is 6.
This indicates the computational complexity without any optimization.

  * •

Optimized scaling: The scaling of the calculation using an optimized approach
is also 6. This indicates the computational complexity after optimization,
showing no reduction in this case.

  * •

Naive FLOP count: The number of floating-point operations (FLOPs) in a naive
approach is $3.266\times 10^{5}$. This measures the computational effort
required without optimization.

  * •

Optimized FLOP count: The number of FLOPs in an optimized approach is
$1.120\times 10^{5}$. This shows the computational effort required after
optimization, indicating a significant reduction.

  * •

Theoretical speedup: The theoretical speedup achieved by optimization is
approximately 2.917 times. This is the ratio of naive FLOP count to optimized
FLOP count, indicating the efficiency gained.

  * •

Largest intermediate: The largest intermediate tensor has $7.776\times 10^{3}$
elements. This represents the size of the largest tensor generated during the
contraction process.

####  VI.6.2 After Decomposition

This time, we performed TT decomposition on the TSP problem. Initially, the
order-6 tensor of size $(6,6,6,6,6,6)$ was decomposed into six tensors with
dimensions $[(6,2),(2,6,3),(3,6,4),(4,6,4),(4,6,2),(2,6)]$. Now, we will use
these tensors to perform contraction with the vector $x$.

This time, we performed TT decomposition on the TSP problem. Initially, the
order-6 tensor of size $(6,6,6,6,6,6)$ was decomposed into six tensors with
dimensions $[(6,2),(2,6,3),(3,6,4),(4,6,4),(4,6,2),(2,6)]$. Now, we will use
these tensors to perform contraction with the vector $x$.

The result of the numpy einsum path after performing the TT decomposition and
contracting the resulting tensors with the vector $x$ is as follows:

    Complete contraction:  i,j,k,l,m,n,iA,AjB,
                           BkC,ClD,DmE,En->
           Naive scaling:  11
       Optimized scaling:  3
        Naive FLOP count:  1.075e+08
    Optimized FLOP count:  7.010e+02
     Theoretical speedup:  153345.826
    Largest intermediate:  1.600e+01 elements

As a result, the maximum scaling parameter is limited to 3, significantly
reduced from 6 before the decomposition. When performing TT decomposition, we
did not apply any approximations and only optimized the computation path. As
we can see, the computation has become significantly lighter after performing
TT decomposition.

##  VII Conclusion

In this study, we introduced HOBOTAN, a Higher Order Binary Optimization
(HOBO) solver that leverages tensor networks and PyTorch to efficiently solve
combinatorial optimization problems. Throughout the paper, we demonstrated the
capabilities of HOBOTAN by solving various complex problems, including the
seating arrangement problem, the search for Pythagorean triples, and the
Traveling Salesman Problem (TSP) using integer encoding.

By utilizing integer encoding, we significantly reduced the number of required
qubits compared to traditional binary encoding methods, such as one-hot
encoding. This reduction in qubits allowed for more efficient computation and
optimized resource usage. Specifically, in the context of the TSP, the integer
encoding of city visit sequences proved to be a powerful approach. This method
enabled us to solve the TSP more efficiently by representing paths and
distances in a compact form, thereby minimizing the computational overhead.

The incorporation of tensor networks into the HOBOTAN framework allowed us to
efficiently manage and optimize high-dimensional tensors. By carefully
designing contraction paths and employing decomposition techniques such as
Singular Value Decomposition (SVD) and Tensor Train Decomposition, we were
able to significantly enhance the computational efficiency of the solver.
These optimizations are critical in handling large-scale problems that would
otherwise be infeasible with conventional methods.

Furthermore, the integration of PyTorch provided powerful GPU acceleration
capabilities, which are essential for handling the parallel computations
involved in tensor contractions. The tensor batch system, introduced in
HOBOTAN, further enhanced the solver’s performance by enabling simultaneous
computation of multiple vector contractions. This system proved particularly
effective in real-world scenarios, allowing HOBOTAN to solve up to 10,000
samples rapidly and efficiently.

In conclusion, HOBOTAN demonstrates a significant advancement in solving
higher-order binary optimization problems by leveraging the combined strengths
of tensor networks and PyTorch. The solver’s ability to efficiently handle
complex problems with reduced computational resources makes it a promising
tool for a wide range of applications in combinatorial optimization. Future
work will focus on expanding HOBOTAN’s capabilities, including multi-GPU
support and further optimization of tensor network techniques, to tackle even
larger and more complex problems.

##  VIII Discussion

The current implementation of the HOBO solver has proven effective in solving
complex optimization problems, confirming its utility. Given that it is based
on tensor networks, a technology not yet fully exploited in this field,
further exploration of this technology is necessary to enhance both speed and
utility. Additionally, the continuous advancements in high-performance
hardware offer a promising outlook for leveraging problem size and speed.

While the basic implementation of HOBOTAN includes numerous optimizations,
there remains room for future improvements, such as supporting multi-GPU
setups to further enhance computational capabilities. This opens up avenues
for scaling the solver to tackle even larger and more complex problems
efficiently.

## References

  * Domino  _et al._ [2022] K. Domino, A. Kundu, . Salehi, and K. Krawiec, Quadratic and higher-order unconstrained binary optimization of railway rescheduling for quantum computing, Quantum Information Processing 21, 10.1007/s11128-022-03670-y (2022).
  * Suksmono and Minato [2022] A. B. Suksmono and Y. Minato, Quantum computing formulation of some classical hadamard matrix searching methods and its implementation on a quantum computer, Scientific Reports 12, 197 (2022).
  * Minato [2024] Y. Minato, Tensor network based hobo solver (2024), arXiv:2407.16106 [quant-ph] .
  * Wang  _et al._ [2023] M. Wang, Y. Pan, Z. Xu, X. Yang, G. Li, and A. Cichocki, Tensor networks meet neural networks: A survey and future perspectives (2023), arXiv:2302.09019 [cs.LG] .
  * Lucas [2014] A. Lucas, Ising formulations of many np problems, arXiv preprint arXiv:1302.5843 (2014).
  * Glos  _et al._ [2022] A. Glos, A. Krawiec, and Z. Zimborás, Space-efficient binary optimization for variational quantum computing, npj Quantum Information 8, 39 (2022).


# Asymmetric Fraunhofer pattern in Josephson junctions from heterodimensional
superlattice V5S8

Juewen Fan(gbsn范珏雯)  Bingyan Jiang(gbsn江丙炎)  Jiaji Zhao(gbsn赵嘉佶)  Ran
Bi(gbsn毕然)  State Key Laboratory for Artificial Microstructure and Mesoscopic
Physics, Frontiers Science Center for Nano-optoelectronics, Peking University,
Beijing 100871, China  Jiadong Zhou(gbsn周家东)  Key Laboratory of Advanced
Optoelectronic Quantum Architecture and Measurement (Ministry of Education),
Beijing Key Laboratory of Nanophotonics & Ultrafine Optoelectronic Systems,
and School of Physics, Beijing Institute of Technology, Beijing 100081, China
Zheng Liu(gbsn刘政)  School of Materials Science and Engineering, Nanyang
Technological University, Singapore 639798, Singapore  Guang Yang(gbsn杨光)  Jie
Shen(gbsn沈洁)  Fanming Qu(gbsn屈凡明)  Li Lu(gbsn吕力)  Beijing National Laboratory
for Condensed Matter Physics, Institute of Physics, Chinese Academy of
Sciences, Beijing 100190, China  Ning Kang(gbsn康宁)  Key Laboratory for the
Physics and Chemistry of Nanodevices and Department of Electronics, Peking
University, Beijing 100871, China  Xiaosong Wu(gbsn吴孝松<EMAIL_ADDRESS>State
Key Laboratory for Artificial Microstructure and Mesoscopic Physics, Frontiers
Science Center for Nano-optoelectronics, Peking University, Beijing 100871,
China  Collaborative Innovation Center of Quantum Matter, Beijing 100871,
China  Shenzhen Institute for Quantum Science and Engineering, Southern
University of Science and Technology, Shenzhen 518055, China

###### Abstract

Introduction of spin–orbit coupling (SOC) in a Josephson junction (JJ) gives
rise to unusual Josephson effects. We investigate JJs based on a newly
discovered heterodimensional superlattice V5S8 with a special form of SOC. The
unique homointerface of our JJs enables elimination of extrinsic effects due
to interfaces and disorder. We observe asymmetric Fraunhofer patterns with
respect to both the perpendicular magnetic field and the current. The
asymmetry is influenced by an in-plane magnetic field. Analysis of the pattern
points to a nontrivial spatial distribution of the Josephson current that is
intrinsic to the SOC in V5S8.

Fraunhofer pattern, Josephson junction, spin–orbit coupling, homointerface

###### pacs:

74.50.+r, 74.45.+c, 03.75.Lm, 71.70.Ej

In the past decade, the interplay between spin–orbit coupling (SOC) and
superconductivity through the proximity effect has received significant
attentionBuzdin2008 ; Reynoso2008Sep ; Lutchyn2010Aug ; Bergeret2015 ; Wu2019
. In materials with strong SOC, a variety of interesting effects may occur
when combining them with superconductors, depending on the form of the SOC,
structural symmetry, etcBuzdin2008 ; Liu2010 ; Yokoyama2014May ;
Yokoyama2014nov ; Bergeret2015 ; Konschelle2015 ; Rasmussen2016 . For
instance, transition metal dichalcogenides with an Ising SOC can turn a
conventional s-wave superconductor into an Ising superconductor through the
proximity effectWu2019 . Another example, which has attracted enormous
interest, is the construction of the Majorana zero mode in spin–orbit coupled
nanowires or topological insulators in close proximity to s-wave
superconductorsLutchyn2010Aug ; Oreg2010Oct ; Fu2008 ; Lyu2018Oct ;
Banerjee2018Dec ; Zheng2019Jun . In Josephson junctions (JJs), the combination
of a Zeeman field and a Rashba spin–orbit interaction can give rise to a
supercurrent, even if the phase difference between two superconductors is
zero, realizing a $\varphi_{0}$ junctionBuzdin2008 ; Reynoso2008Sep ;
Bergeret2015 .

A major playground for such studies is superconductor–normal
metal–superconductor (SNS) JJs, in which the normal metal provides SOC. It is
well-known that the critical current of a conventional JJ as a function of
magnetic field exhibits a Fraunhofer interference pattern, described
byTinkham1996

$I_{\mathrm{c}}(B_{z})=I_{\mathrm{c}}(0)\left|\frac{\sin(\pi
B_{z}S/\Phi_{0})}{\pi B_{z}S/\Phi_{0}}\right|$ (1)

, where $I_{\mathrm{c}}$ is the critical current, $B_{z}$ the magnetic field
perpendicular to the junction plane, $S$ the effective area of junction, and
$\Phi_{0}=h/2e$ is the flux quantum. $I_{\mathrm{c}}$ is apparently
independent of the polarity of $B_{z}$,
$I_{\mathrm{c}}^{\pm}(B_{z})=I_{\mathrm{c}}^{\pm}(-B_{z})$, where $\pm$
indicates the current polarity. Its amplitude usually remains unchanged upon
reversing the current $I$, too,
$I_{\mathrm{c}}^{+}(B_{z})=-I_{\mathrm{c}}^{-}(B_{z})$. That is, the
Fraunhofer pattern is symmetric in $B_{z}$ and $I$. However, recent studies
have shown that these symmetries can be broken, especially in the presence of
SOC. Broken current symmetry may be induced by structural asymmetry of the
junctionHu2007Aug ; Wu2021Mar ; Misaki2021Jun , or a multiband
effectReynoso2008Sep ; Yokoyama2014May ; Trimble2021Jun . Broken field
symmetry was attributed to disorder potentialsSuominen2017 , or nonuniform
SOCChen2018a ; Assouline2019 ; Beach2021 . Most of these cases require
extrinsic mechanisms, which often exist in a practical junction. A general
relation between the pattern asymmetry and the combination of a spin–orbit
coupling, in-plane magnetic field and disorder has been theoretically
studiedLiu2010 ; Rasmussen2016 .

Here we study JJs fabricated with a novel material,
$\mathrm{V}_{5}\mathrm{S}_{8}$ of a heterodimensional superlattice, whose
special form of SOC gives rise to a rare in-plane Hall effectZhou2021 . The
homointerface of the junction enables investigation of intrinsic effects. It
is found that both the magnetic field and current symmetries of the Fraunhofer
pattern are broken. Moreover, the pattern can be tuned by an in-plane field.
Consistent patterns and field dependence are observed in all devices. The
phenomenon suggests a nontrivial spatial distribution of the Josephson current
that is intrinsic to V5S8 and stems from its SOC.

Micro strips of V5S8 films were grown by chemical vapor deposition on SiO2
substratesZhou2021 . Strips with about 10 nm thickness were selected for this
study. As for superconductors, 31 nm aluminum, following a 2 nm titanium
adhesion layer, were e-beam deposited. Some measurements were performed in an
Oxford 3He cryostat with a base temperature of 250 mK. Experiments at a lower
temperature were carried out in an Oxford dilution refrigerator with a base
temperature of 10 mK. The magnetic field was provided by a three-axis vector
superconducting magnet. To avoid the effect of hysteresis due to heating
during a bias current sweep, the sweep always started from zero.

V5S8 film used in this study has a unique superlattice structure, in sharp
contrast with conventional V5S8 compoundsZhou2021 . It can be seen as VS2
layers intercalated with V2S2 atomic chains. It exhibits an extraordinary in-
plane Hall effect, suggesting a peculiar form of spin–orbit-coupling. Because
Al(31 nm)/V5S8(10 nm) bilayers are non-superconducting due to a strong inverse
proximity effect, aluminum JJs can be made using the Al(31 nm)/V5S8(10 nm)
bilayer as the normal metal weak linkFan2021 . The JJ is a cross-bar structure
made from a V5S8 strip and an aluminum strip. The resultant homointerface
between the superconductor and normal metal is highly transparent and uniform,
providing an ideal JJ platform for investigation of the intrinsic effect of
SOC.

The measured JJ device is illustrated in Fig. 1a. Since our superlattice V5S8
is a polar crystal, a coordinate system can be uniquely defined, which is
necessary for the latter discussion. Note that the growth of V5S8 strips
starts from a nucleation center and propagates along the $y$ axis. It is thus
reasonable to assume that the growth is consistently towards the same
crystalline direction. This direction is defined as the $-y$ direction. A
positive current flows in the $+x$ direction of the crystal structure, while a
positive magnetic field is in the $+z$ direction, which is pointing out of the
plane of the substrate. The differential resistance $\mathrm{d}V/\mathrm{d}I$
displays a Fraunhofer-like pattern in the perpendicular magnetic field $B_{z}$
versus bias current $I$ mapping at 250 mK (Fig. 1b). The period of the pattern
$B_{z}^{\mathrm{p}}$ is 0.99 mT, yielding an effective junction area
$S=\Phi_{0}/B_{z}^{\mathrm{p}}=2.09$ $\mu\mathrm{m}^{2}$, which is close to
the nominal area of the JJ $S_{\mathrm{n}}=1.61$ $\mu\mathrm{m}^{2}$. The
discrepancy comes from the London penetration depth and flux focusing. The
diffusion constant of junction, $D$, is 6.9$\times 10^{-3}$ m${}^{2}\cdot$s-1,
calculated by $D=\frac{1}{3}(\frac{\pi
k_{\mathrm{B}}}{e})^{2}\frac{\sigma}{\gamma}$Pippard1960 . Here,
$\sigma=3.9\times 10^{7}$ S$\cdot$m-1 is the electrical conductivity of the
weak link, $\gamma=1.4\times 10^{2}$ J$\cdot$m${}^{-3}\cdot$K-2 is the
electronic specific heat coefficient of AlKittel2005 . Since the electrical
conductivity of Al film is much larger than that of V5S8, the electronic
specific heat coefficient of the weak link can be approximate to that of Al.
Then the mean free path $l_{\mathrm{e}}$ is 10.2 nm through
$D=v_{\mathrm{F}}l_{\mathrm{e}}/3$, where $v_{\mathrm{F}}=2.02\times 10^{8}$
cm$\cdot$s-1 is the Fermi velocity of AlKittel2005 . Because $l_{\mathrm{e}}$
is much shorter than the junction length $L=0.92$ $\mu$m, the junction is in
the diffusive limit, in which case the Thouless energy $E_{\mathrm{Th}}$ is
given by $\hbar D/L^{2}$Dubos2001  and becomes 5.4 $\mu$eV, much smaller than
the superconducting gap, indicating that the JJ is in the long junction limit.

Figure 1: Fraunhofer interference patterns of a JJ. (a) An optical micrograph
of the device with an illustration of the measurement configuration. The
junction length is $L=0.92$ $\mu$m, the junction width is $W=1.75$ $\mu$m. The
scale bar is 1 $\mu$m. (b) The interference pattern at 250 mK. The dashed
white lines indicate the position of the nodal points and $I=0$. (c) The zoom-
up pattern taken at 32 mK. (d) The further zoom-up scan on the first-order
nodes in (c). A break in the $x$ axis is introduced for better comparison.

The critical current at the low-order nodal points does not vanish as
anticipated by a standard Fraunhofer pattern, which will be discussed later.
The most unusual feature is that the pattern is lack of the mirror symmetries
with respect to both the current and the field,
$I_{\mathrm{c}}^{+}(B_{z})\neq-I_{\mathrm{c}}^{-}(B_{z})$ and
$I_{\mathrm{c}}^{\pm}(B_{z})\neq I_{\mathrm{c}}^{\pm}(-B_{z})$. The nodal
points in the upper half-plane ($I^{+}$) are not aligned with those in the
lower half-plane ($I^{-}$). The misalignment can be seen more clearly in a
zoom up mapping taken at low temperatures, where the superconducting
transition is sharp. As depicted in Fig. 1c, the nodal points in the upper
half-plane shift to the right with respect to those in the lower half-plane. A
further zoom up scan on the first-order nodes shows that aside from the
horizontal shift, the nodal points also shift vertically (Fig. 1d). These in
the second and forth quadrants, assigned numbers of $n=2$ and 4, respectively,
are closer to the $I=0$ axis than those in the first and third quadrants,
assigned numbers of $n=1$ and 3. In all of 7 devices we have measured, the
nodal points shift towards the same direction, suggesting that the effect is
intrinsic to V5S8, rather than due to the sample geometry, differences at the
interface or disorder. On the other hand, the inversion symmetry
$-I_{\mathrm{c}}^{+}(B_{z})=I_{\mathrm{c}}^{-}(-B_{z})$ is retained.

We now apply an in-plane field and investigate the change in the nodal point.
Fig. 2 shows 12 fine mappings around the nodal points in Fig. 1d in the
presence of the in-plane field of $B_{\parallel}=0.12$ T along different
directions. Let $\theta$ denote the angle between $B_{\parallel}$ and the $x$
axis. One may immediately notice that the inversion symmetry of pattern is now
broken. The critical currents at four nodal points, $I_{\mathrm{c}}^{n}$
$(n=1,2,3,4)$, are plotted as a function of the field angle $\theta$, shown in
Fig. 3a. $I_{\mathrm{c}}^{n}$s oscillate with a $2\pi$ period.
$|I_{\mathrm{c}}^{1}(\theta)|$ and $|I_{\mathrm{c}}^{3}(\theta)|$ are exactly
out-of-phase, so are $|I_{\mathrm{c}}^{2}(\theta)|$ and
$|I_{\mathrm{c}}^{4}(\theta)|$. Only at two critical angles,
$\theta=158^{\circ}$ and $338^{\circ}$, $|I_{\mathrm{c}}^{1}|$ is equal to
$|I_{\mathrm{c}}^{3}|$ and $|I_{\mathrm{c}}^{2}|$ is equal to
$|I_{\mathrm{c}}^{4}|$, indicating restoration of the inversion symmetry. To
quantify the degree of inversion symmetry breaking, let us employ the
difference of $I_{\mathrm{c}}$ between a pair of nodes that are inversion
symmetric in the absence of an in-plane field. As illustrated in Fig. 3b, the
asymmetry is at its maximum at $\theta=68^{\circ}$ and $248^{\circ}$.
Similarly, we can assess the asymmetry from the field of nodal points. Let
$B_{z}^{n}$ $(n=1,2,3,4)$ denote the perpendicular field at the nodal point.
To eliminate the estimation error of the zero field due to unavoidable sample
misalignment, we calculate $\Delta B_{z}^{23}=B_{z}^{2}-B_{z}^{3}$ and $\Delta
B_{z}^{14}=B_{z}^{1}-B_{z}^{4}$, plotted as a function of $\theta$ in Fig. 3c.
The more $\Delta B_{z}^{23}$ deviates from $\Delta B_{z}^{14}$, the more
asymmetric the pattern is. The evolution of the asymmetry is consistent with
that inferred from $I_{\mathrm{c}}$. Comprehensive measurements of the angular
dependence of the pattern on another sample, shown in the supplementary
material, yield an almost identical behavior and a similar critical angle at
which $|I_{\mathrm{c}}^{1}|$($|I_{\mathrm{c}}^{2}|$) intercepts
$|I_{\mathrm{c}}^{3}|$($|I_{\mathrm{c}}^{4}|$).

Figure 2: Fine mapping of the first-order nodes in the interference pattern in
the presence of $B_{\parallel}=0.12$ T oriented at different angles with
respect to the $x$ axis. Two dashed yellow lines at $I=\pm 3.158$ $\mu$A are
used as a baseline to highlight the shift of the critical current at the nodal
points. Figure 3: Angular dependence of the position of the nodal points. (a)
Absolute value of the critical current. (b) The difference of the critical
currents, $|I_{\mathrm{c}}^{1}|-|I_{\mathrm{c}}^{3}|$ and
$|I_{\mathrm{c}}^{2}|-|I_{\mathrm{c}}^{4}|$. The dashed curves represent fits
to a sine function. (c) The difference of the field at the nodal points,
$B_{z}^{2}-B_{z}^{3}$ and $B_{z}^{1}-B_{z}^{4}$.

Since $|I_{\mathrm{c}}^{1}(\theta)|$ and $|I_{\mathrm{c}}^{3}(\theta)|$
oscillate with the same amplitude and opposite phase, i.e.,
$|I_{\mathrm{c}}^{1}(\theta)|=|I_{\mathrm{c}}^{3}(\theta+\pi)|$, we have
$I_{\mathrm{c}}^{1}(B_{\parallel})=-I_{\mathrm{c}}^{3}(-B_{\parallel})$,
which, combined with the similar behavior of $\Delta B_{z}^{23}$ and $\Delta
B_{z}^{14}$, suggests

$I_{\mathrm{c}}^{+}(B_{z},B_{\parallel})=-I_{\mathrm{c}}^{-}(-B_{z},-B_{\parallel})$
(2)

. In another word, the diffraction pattern is inversion symmetric if the total
magnetic field is considered. This can also be directly seen in Fig. 2, where
the pattern at $\theta$ is inversion symmetric to that at
$\theta+180^{\circ}$.

First, let us discuss the reason for the symmetry breaking at zero in-plane
field seen in Fig. 1. One possible cause is the self-field effectBarone1982 .
For JJs with certain geometries, a nonuniform current distribution in the
superconductor electrodes can produce an additional field through the junction
area, which shifts the pattern horizontally. As the self-field depends on the
polarity of the current, the patterns in the upper and lower half-plane shift
towards opposite directions. The geometry of our JJs has a mirror symmetry
about the center line along the current direction, precluding the self-field
effect. Even if the device slightly deviates from the mirror symmetry due to
imperfections introduced by fabrication, the skewness of the pattern should
vary from device to device. The consistent shift in all devices indicates that
the symmetry breaking is associated with the crystal structure of V5S8.

Abrikosov vortices in superconducting electrodes near the JJ can lead to an
asymmetric Fraunhofer pattern, too. However, the effect breaks the field
symmetry, $I_{\mathrm{c}}^{\pm}(B_{z})=I_{\mathrm{c}}^{\pm}(-B_{z})$, leaving
the current symmetry, $I_{\mathrm{c}}^{+}(B_{z})=-I_{\mathrm{c}}^{-}(B_{z})$,
intactGolod2010 ; Golod2019 , which is at odds with our results. Moreover, we
intentionally made the superconductor leads narrow enough so that vortices can
not enterStan2004 . For comparison, when the leads are wider than 0.65 $\mu$m,
jittering of the pattern due to vortex jump started to appear in the
Fraunhofer pattern.

Rasmussen et al. have studied the interference pattern of a SNS JJ
theoretically and related its symmetry to certain symmetries in the
HamiltonianRasmussen2016 . The symmetries that are required to break in order
to break the mirror symmetries,
$I_{\mathrm{c}}^{\pm}(B_{z})=I_{\mathrm{c}}^{\pm}(-B_{z})$ and
$I_{\mathrm{c}}^{+}(B_{z})=-I_{\mathrm{c}}^{-}(B_{z})$, are listed in the
Table. 1. Here $\sigma_{x,y,z}$ are the Pauli matrices, $P_{x(y)}$ is the
$x(y)$-parity operation, $T$ is the time reversal operation. Disordered
potential $V_{x(y)}$, which is asymmetric under $P_{x(y)}$, unlikely plays an
important role, otherwise the pattern in different devices would shift in
random directionsSuominen2017 , not to mention the fact that our JJs, made
from a continuous aluminum film, are presumably quite uniform. In the absence
of an in-plane field, two forms of SOC, Rashba SOC $\alpha$ and Dresselhaus
SOC $\beta$, can break most of the symmetries. It is plausible that V5S8 has
both $\alpha$ and $\beta$ or similar forms of SOC, as those SOCs are the key
ingredients in the theories on the in-plane Hall effectZhang2011 ; Ren2016a ;
Liu2018a ; Zyuzin2020 , which has been observed in V5S8. However, there are
still two more symmetry operators left, $\sigma_{z}P_{y}P_{x}$ and
$\sigma_{z}P_{x}P_{y}T$, which guarantee the mirror symmetries of the pattern.

Table 1: Symmetry operations $U$ protecting $-I_{\mathrm{c}}^{+}(B_{z})=I_{\mathrm{c}}^{-}(B_{z})$ and $I_{\mathrm{c}}^{\pm}(B_{z})=I_{\mathrm{c}}^{\pm}(-B_{z})$ respectively adapted with permission from Ref. Rasmussen2016, .Copyrighted by the American Physical Society. $-I_{\mathrm{c}}^{+}(B_{z})=I_{\mathrm{c}}^{-}(B_{z})$ | $I_{\mathrm{c}}^{\pm}(B_{z})=I_{\mathrm{c}}^{\pm}(-B_{z})$
---|---
$U$ |  Broken by |  $U$ |  Broken by
$P_{y}P_{x}$ |  $\alpha$, $\beta$, $V_{x}$, $V_{y}$ |  $\sigma_{x}P_{y}$ |  $B_{y}$, $\alpha$, $V_{y}$
$\sigma_{z}P_{y}P_{x}$ |  $B_{x}$, $B_{y}$, $V_{x}$, $V_{y}$ |  $\sigma_{y}P_{y}$ |  $B_{x}$, $\beta$, $V_{y}$
$\sigma_{x}P_{y}T$ |  $B_{x}$, $\alpha$, $V_{y}$ |  $P_{x}P_{y}T$ |  $B_{x}$, $B_{y}$, $\alpha$, $\beta$, $V_{x}$, $V_{y}$
$\sigma_{y}P_{y}T$ |  $B_{y}$, $\beta$, $V_{y}$ |  $\sigma_{z}P_{x}P_{y}T$ |  $V_{x}$, $V_{y}$

It has been proposed that the mirror symmetries of the interference pattern,
$I_{\mathrm{c}}^{\pm}(B_{z})=I_{\mathrm{c}}^{\pm}(-B_{z})$ and
$I_{\mathrm{c}}^{+}(B_{z})=-I_{\mathrm{c}}^{-}(B_{z})$, are broken in SNS
based on inversion symmetry breaking topological materialsChen2018 . This is
because the topological surface states on two sides of the JJ have different
Fermi velocities, resulting in different Josephson currents. The interference
pattern by the surface states is thus shifted and skewed, depending on the
polarity of the current. Taken into account the dominant contribution from the
bulk, the resultant inversion symmetric pattern agrees with what we observed.
Broken mirror symmetries of the Fraunhofer pattern were observed in
topological materialsBocquillon2016 ; Kononov2020Jun . The difference of the
Josephson currents on two sides breaks the mirror symmetry of the Josephson
current spatial distribution about the center line of JJ along the current. It
is worth noting that this mirror symmetry is essential for the nodal point
being on the $I=0$ axis. Therefore, the finite critical current of the nodal
point in our data provides additional evidence, independent of the model on
inversion symmetry breaking topological materials, that the Josephson currents
on the two sides are not equal.

The asymmetry of the pattern is tuned by the in-plane field, seen in Fig. 2.
It has been shown that an in-plane field can suppress the Josephson current in
the center of a JJ via a peculiar flux focusing effectSuominen2017 .
Consequently, the zero field critical current is reduced and the interference
pattern becomes more like that of a superconducting quantum interference
device. It works when the field is parallel to the current and does not depend
on the field polarity. As the field is rotating in the plane, a period of
$\pi$ is expected in the angular dependence of the critical current, in
contrast to a period of $2\pi$ shown in Fig. 3. Furthermore, the effect does
not depend on the polarity of the current, so it will not break aforementioned
two mirror symmetries of the pattern. At last, the flux focusing effect
maximizes at $\theta=0^{\circ}$ and $180^{\circ}$, whereas the critical
currents of the nodal points at these two angles are at neither the maximum
nor the minimum, shown in Fig. 3. These distinctions rule out the flux
focusing effect of the in-plane field.

An in-plane field usually acts on electrons by the Zeeman effect, thus its
strong influence on the pattern suggests that the asymmetry is related to the
spin. In the case of inversion symmetry broken topological materials, the
Zeeman effect will break the time reversal symmetry of the surface states,
which alters the Josephson current on the two sides. The change in the current
distribution inevitably leads to a change of the pattern. The critical angles
of $\theta=158^{\circ}$ and $338^{\circ}$, at which the inversion symmetry of
the pattern is restored, are probably related to the particular form of SOC in
V5S8.

The unusual symmetry in the Fraunhofer pattern of Al–(Al/V5S8)–Al Josephson
junctions suggests an anomalous Josephson coupling stemming from the SOC in
heterodimensional superlattice V5S8. This newly discovered compound, with its
exotic SOC and the ability of forming a homointerface Josephson junction,
offers an intriguing platform for investigation of unconventional Josephson
effects and superconductivity.

## Supporting Information

The measurements of the in-plane field dependence of the Fraunhofer pattern on
another Al–(Al/V5S8)–Al Josephson junction S2.

###### Acknowledgements.

Project supported by the National Key Basic Research Program of China (Grant
No. 2016YFA0300600) and the National Natural Science Foundation of China
(Grant Nos. 11574005 and 11774009).

## References

  * (1) Buzdin A 2008 Phys. Rev. Lett. 101 107005
  * (2) Reynoso A A, Usaj G, Balseiro C A, Feinberg D and Avignon M 2008 Phys. Rev. Lett. 101 107001
  * (3) Lutchyn R M, Sau J D and Das Sarma S 2010 Phys. Rev. Lett. 105 077001
  * (4) Bergeret F S and Tokatly I V 2015 Europhys. Lett. 110 57005
  * (5) Wu Y, He J J, Han T, Xu S, Wu Z, Lin J, Zhang T, He Y and Wang N 2019 Phys. Rev. B 99 121406
  * (6) Liu J F and Chan K S 2010 Phys. Rev. B 82 125305
  * (7) Yokoyama T, Eto M and Nazarov Y V 2014 Phys. Rev. B 89 195407
  * (8) Yokoyama T and Nazarov Y V 2014 Europhys. Lett. 108 47009
  * (9) Konschelle F, Tokatly I V and Bergeret F S 2015 Phys. Rev. B 92 125443
  * (10) Rasmussen A, Danon J, Suominen H, Nichele F, Kjaergaard M and Flensberg K 2016 Phys. Rev. B 93 155406
  * (11) Oreg Y, Refael G and von Oppen F 2010 Phys. Rev. Lett. 105 177002
  * (12) Fu L and Kane C L 2008 Phys. Rev. Lett. 100 096407
  * (13) Lyu Z, Pang Y, Wang J, Yang G, Fan J, Liu G, Ji Z, Jing X, Yang C, Qu F and Lu L 2018 Phys. Rev. B 98 155403
  * (14) Banerjee A, Sundaresh A, Ganesan R and Kumar P S A 2018 ACS Nano 12 12665
  * (15) Zheng H and Jia J F 2019 Chin. Phys. B 28 67403
  * (16) Tinkham M 1996 Introduction to Superconductivity (New York: McGraw-Hill) pp. 215–217
  * (17) Hu J, Wu C and Dai X 2007 Phys. Rev. Lett. 99 067004
  * (18) Wu H, Wang Y, Xu Y, Sivakumar P K, Pasco C, Filippozzi U, Parkin S S P, Zeng Y J, McQueen T and Ali M N 2021 arXiv:2103.15809 [cond-mat.supr-con]
  * (19) Misaki K and Nagaosa N 2021 Phys. Rev. B 103 245302
  * (20) Trimble C J, Wei M T, Yuan N F Q, Kalantre S S, Liu P, Han H J, Han M G, Zhu Y, Cha J J, Fu L and Williams J R 2021 npj Quantum Mater. 6 61
  * (21) Suominen H J, Danon J, Kjaergaard M, Flensberg K, Shabani J, Palmstrøm C J, Nichele F and Marcus C M 2017 Phys. Rev. B 95 035307
  * (22) Chen A Q, Park M J, Gill S T, Xiao Y, Reig-i Plessis D, MacDougall G J, Gilbert M J and Mason N 2018 Nat. Commun. 9 3478
  * (23) Assouline A, Feuillet-Palma C, Bergeal N, Zhang T, Mottaghizadeh A, Zimmers A, Lhuillier E, Eddrie M, Atkinson P, Aprili M and Aubin H 2019 Nat. Commun. 10 126
  * (24) Beach A, i Plessis D R, MacDougall G and Mason N 2021 J. Phys.: Condens. Matter 33 425601
  * (25) Zhou J, Zhang W, Lin Y C, Zhou Y, Du H, Tang B, Shi J, Jian B, Cao X, Lin B, Zhu C, Deng Y, Fu Q, Duan R, Wang X, Chen J, Guo S, Guo W, Huang Y, Yao Y, Gao Y, Yao Y, Suenaga K, Wu X S and Liu Z 2021 Heterodimensional superlattice with room temperature anomalous Hall effect, under review
  * (26) Fan J, Jiang B, Zhao J, Bi R, Zhou J, Liu Z, Kang N, Qu F, Lu L and Wu X 2022 J. Appl. Phys. 131 093903
  * (27) Pippard A B 1960 Rep. Prog. Phys. 23 176
  * (28) Kittel C 2005 Introduction to Solid State Physics (New Jersey: John Wiley & Sons) pp. 146,139
  * (29) Dubos P, Courtois H, Pannetier B, Wilhelm F K, Zaikin A D and Schön G 2001 Phys. Rev. B 63 064502
  * (30) Barone A and Paternò G 1982 Physics and Applications of the Josephson Effect (New York: John Wiley & Sons) pp. 96–112
  * (31) Golod T, Rydh A and Krasnov V M 2010 Phys. Rev. Lett. 104 227003
  * (32) Golod T, Pagliero A and Krasnov V M 2019 Phys. Rev. B 100 174511
  * (33) Stan G, Field S B and Martinis J M 2004 Phys. Rev. Lett. 92 097003
  * (34) Zhang Y and Zhang C 2011 Phys. Rev. B 84 085123
  * (35) Ren Y, Zeng J, Deng X, Yang F, Pan H and Qiao Z 2016 Phys. Rev. B 94 085411
  * (36) Liu Z, Zhao G, Liu B, Wang Z F, Yang J and Liu F 2018 Phys. Rev. Lett. 121 246401
  * (37) Zyuzin V A 2020 Phys. Rev. B 102 241105
  * (38) Chen C Z, He J J, Ali M N, Lee G H, Fong K C and Law K T 2018 Phys. Rev. B 98 075430
  * (39) Bocquillon E, Deacon R S, Wiedenmann J, Leubner P, Klapwijk T M, Brüne C, Ishibashi K, Buhmann H and Molenkamp L W 2016 Nat. Nanotechnol. 12 137
  * (40) Kononov A, Abulizi G, Qu K, Yan J, Mandrus D, Watanabe K, Taniguchi T and Schönenberger C 2020 Nano Lett. 20 4228

Supplemental Materials

This Supplemental Material Section contains the measurements of the in-plane
field dependence of the Fraunhofer pattern on another Al–(Al/V5S8)–Al
Josephson junction S2.

Figure S1: Fine mapping of the first-order nodes in the interference pattern
of JJ S2 in the presence of $B_{\parallel}=0.12$ T oriented at different
angles with respect to the $x$ axis. Two dashed yellow lines at $I=\pm 0.9$
$\mu$A are used as a baseline to highlight the shift of the critical current
at the nodal points.

JJ S2 in Fig. S2 displays a behavior similar to that of the JJ in the main
text, including the critical angles at which
$|I_{\mathrm{c}}^{1}|$($|I_{\mathrm{c}}^{2}|$) intercepts
$|I_{\mathrm{c}}^{3}|$($|I_{\mathrm{c}}^{4}|$). For JJ S2, the critical angles
are about $140^{\circ}$ and $320^{\circ}$.

Figure S2: Angular dependence of the position of the nodal points of JJ S2.
(a) Absolute value of the critical current. (b) The difference of the critical
currents, $|I_{\mathrm{c}}^{1}|-|I_{\mathrm{c}}^{3}|$ and
$|I_{\mathrm{c}}^{2}|-|I_{\mathrm{c}}^{4}|$. The dashed curves represent fits
to a sine function. (c) The difference of the field at the nodal points,
$B_{z}^{2}-B_{z}^{3}$ and $B_{z}^{1}-B_{z}^{4}$.