Long range magnetic dipole-dipole interaction mediated by a superconductor

Quantum computation and simulation requires strong coherent coupling between qubits, which may be spatially separated. Achieving this coupling for solid-state based spin qubits is a long-standing challenge. Here we theoretically investigate a method for achieving such coupling, based on superconducting nano-structures designed to channel the magnetic flux created by the qubits. We detail semi-classical analytical calculations and simulations of the magnetic field created by a magnetic dipole, depicting the spin qubit, positioned directly below nanofabricated apertures in a superconducting layer. We show that such structures could channel the magnetic flux, enhancing the dipole-dipole interaction between spin qubits and changing its scaling with distance, thus potentially paving the way for controllably engineering an interacting spin system.


INTRODUCTION
Solid-state qubits have emerged as a potential quantum information processing architecture, with leading candidates such as atomic defects in bulk materials [1] and quantum dots [2]. High fidelity quantum control of such individual spin qubits has been demonstrated in various systems [3,4]. However, scalable coherent coupling of these qubits, required in order to produce robust twoqubit gates, still poses a significant challenge. The direct magnetic coupling between spin qubits via the dipoledipole interaction is relevant only for spins that are quite close (at the scale of 10 nm), as this interaction usually decays with the distance cubed. This has been demonstrated, e.g., for both Nitrogen-Vacancy (NV) centers in diamond [5,6] and for quantum dots [7,8]. Applications in quantum computing will require a high number of coupled qubits which can be easily addressed. This would be much easier to implement if the qubits are spatially separated to much longer distances than usually possible with dipole-dipole interaction. One such method is to use another system, such as photons, as a quantum bus [2,9,10]. Another method is to use an external cavity or a "floating gate" [11,12].
In this work we propose a method to increase the dipole-dipole interaction to much longer distances. We propose to place dipoles inside fabricated nano-structures in a thin layer of a superconductor. The superconductor would guide the field lines created by the dipoles into the nano-structures, a phenomena known as flux focusing, and increase the interaction between the dipoles. A possible structure would be a "dog-bone" shaped structure shown in Fig. 1. Fig. 1(a) is a side-way view of "dogbone" structure fabricated directly above an NV center or another bulk-semiconductors based qubit. For quantum dots or other dipoles, this device could be fabricated on a Si wafers with the quantum dots later placed inside the apertures. Fig. 1(b) shows a top-down view of a "dogbone" structure. This structure, with a ferromagnet instead of a superconductor, has already been proposed as a possible way to create long-range coherent interaction between NV centers [13]. It has also been proposed as a long-distance coupler between phosphorous based qubits in Si [3]. We also note that this device could be used for sensing application, with the sensor placed in one aperture and the sample placed in another aperture. This could be useful whenever the sample requires a "clean" area far from the sensor. In the case of NV centers, such a realistic scenario could arise for a sample which might be damaged by the laser used to address the NV center.
We look at three simpler scenarios illustrated in Fig. 2: A round aperture with a dipole in the center, a round aperture with a dipole next to the left edge, and an elliptical aperture with a dipole next to the left edge. In all three cases the 2nd dipole would be placed next to the right edge of the aperture. We note that the difference between the "dog-bone" behavior and the ellipse behavior should be small when the ellipse eccentricity is close to 1, due to the similar one-dimensional confinement of the magnetic field. Having a circular aperture at the edge of the ellipse (a "dog-bone" structure) would ease the fabrication and the accurate placement of the structure relative to the qubit, but should not change the power law scaling of the coupling. We show that by placing the dipoles inside these fabricated nano-structures the dipole-dipole interaction is significantly enhanced.
Here we present a method for calculating the exact analytical solution for the magnetic field arising from a single dipole inside a circular aperture in a superconducting thin film, under the assumption of zero penetration depth. We solve for a dipole at the center of the aperture and then at a shifted position, and show that the flux focusing and field confinement can be used to dramatically increase the dipole-dipole interaction. In addition, we numerically solve the problem with a finite penetration depth and show a very good agreement with the analytic calculations. Finally, we extend the numerical results to an elliptical aperture and show an even further increase in the interaction strength.
During the preparation of this article we became aware of a recent similar work used to find the flux focusing in a parallel SQUID array [30].

ANALYTICAL ANALYSIS
A known form of the London Equations, which describes the relation between the vector potential A and the current density J in a superconductor is [31]  where µ 0 is the vacuum permeability and λ is the superconductor penetration depth. Thus λ → 0 ⇒ A| σ → 0, where σ denotes the boundary of the superconductor, implying that the limit of zero penetration depth of a superconductor results in a vanishing vector potential on the boundaries of the superconductor [32]. Under these conditions, we can base our calculation off of the method of the Dirichlet electrostatic Green function [33]. In the Coulomb Gauge we have ∇ 2 A(r) = −µ 0 J (r), where each of the Cartesian components of the vector potential behaves like an electrostatic potential with ρ(r) 0 → µ 0 J i (r) ( 0 is the vacuum permittivity). Using this method, when given a current density J (r) and a Green function G(r, r ), one is able to calculate the vector potential and the resulting magnetic field: In our case, we have a superconducting film lying in the xy plane, with a circular aperture of radius R, i.e.: σ = {z = 0 ∩ ρ ≥ R}. The current density of a point magnetic dipole m, located at r 0 , is given by [33] J (r) = −m × ∇δ(r − r 0 ).
The problem of finding the Green Function of a semiinfinite film, {z = 0 ∩ x ≥ 0}, is treated in [34]. Following the process in [35], one is able to apply the Kelvin inversion transformation that can be used to generate new solutions for the Poisson equation from known ones [36]. This transformation inverts the space relative to a predefined sphere. By applying the Kelvin inversion we are able to transform the Green function of the semi-infinite film to one of an infinite film with a circular aperture, i.e for σ. For more details, see appendix.
Attempting to generalize this approach in order to solve an asymmetric aperture (more specifically an elliptic aperture) is not a trivial task. Naively stretching one of the x or y coordinates in the circular Green function generates a new function which vanishes outside of an elliptical aperture, but no longer satisfies the Poisson equation. There is a modified Kelvin inversion transformation, which inverts the space relative to an ellipsoid rather than a sphere, but unfortunately, this transformation is non-conformal [37]. Non-conformality implies that it will not generate a correct Poisson equation solution in the new elliptical geometry. There are many known 2D conformal transformations that are able to transform a circular aperture to an asymmetric aperture, but applying them to the x, y coordinates of a 3D Green function will not be conformal. To the best of our knowledge there is no 3D conformal transformation which can transform a circle into an ellipse. We therefore proceed with the analysis of a circular aperture, with a centered and off-center dipole, and then continue with numerical calculations for the elliptic case.

Dipole at the center of an aperture
The simplest case to solve is the case illustrated in Fig. 2(a): a dipole in the center of a round aperture. Plugging in the Green function from [35] and the current density of a dipole [33] located at the origin into Eq. 2, we get (For a detailed derivation, see appendix): Here, It is instructive to introduce the vectorn(r), so one can easily notice that the vector potential In Eq. 3 has the same form as a free magnetic dipole vector potential, withn(r) →r. The entire effect of the superconductor is reflected in the coefficient C(r) in theρ direction, which is a location dependent scaling factor. We see that on the boundary σ, α = 0, so n| σ = 0, and then A| σ = 0, consistent with our construction of G(r, r ). As a sanity check, the limit for an infinitely large aperture R → ∞ (or equivalently, reproduces the expression of a free dipole. Now for the magnetic field, taking the curl of Eq. 3 we obtain When considering the limit of an infinitely large aperture, n(r) →r as before, the two right terms cancel each other out and m · ∇ n(r) → m, resulting in the known expression for the magnetic dipole field. In Fig. 3 we plot the magnetic field calculated from Eq. 4, comparing the case of a small aperture (normalized size R = 1) to the free dipole case (without an aperture). The stream plot provides the direction of the field, and the contour plot is scaled relative to the magnitude of the field. We see the comparison to the free magnetic dipole field in the right side plot. It is evident visually that the superconductor confines the magnetic flux.
Restricting ourselves to the xy plane, we can derive a relatively simple analytical expression for the magnetic field inside the aperture. Due to the superconductor fully blocking the field in theẑ direction, the strongest confinement occurs for a dipole oriented along theẑ direction (See appendix for the full derivation). Starting from Eq. 4 and setting the dipole to point along theẑ direction, we derive the magnetic field in the xy plane as This equation is plotted in Fig. 4, stressing the flux focusing and thus enhanced magnetic field near the edge of the aperture. We can now use Eq. 5 to infer the strength of the field confinement near the edge, if we maintain a constant distance from the edge, denoted by d, and increase the aperture's radius: We find that the power-law scaling of the magnetic field with the radius has increased from the usual dipole power-law of −3 to −2.5. This result is plotted in Fig.(5). The field vanishes on the superconductor, since no field lines can cross it. Close to the center, for which the limit of ρ → 0 (or R → ∞) is valid , we find a behavior consistent with a free dipole, approaching − µ 0 m 4πρ 3 . Closer to the edge we see that the field is approaching infinity as can be seen from Eq. 6 Dipole at the side of an aperture Shifting to a dipole which is off-center, e.g. at the left edge of a round aperture [Fig. 2(b)] is relatively simple. Starting with Eq. 2, we use a dipole current density which is shifted in thex direction by x 0 . The derivation is similar to that of the non-shifted case but the solutions for the vector potential and the field, calculated with Mathematica [38], are given as complicated expressions which are too long to include fully. Nevertheless, we verify that by substituting x 0 = 0 we again obtain the results of a centered dipole (Eq. [3][4][5]. Moreover, the magnetic field can be plotted and simplified next to the edge: we derive the expression for a dipole shifted to x 0 = d − R (a distance d from the left edge of the aperture) and we look at the magnetic field symmetrically near the other edge, at the point x = R − d, where d is assumed small relative to the radius R. Using Mathematica [38] we obtain a simplified series expansion We find a significant enhancement of the magnetic field, with the power-law scaling improving from the usual dipole power-law of −3 to −2. This change in scaling is depicted in Fig. 5.
In order to gain insight into this behavior, we can start by examining a dipole which is shifted to x 0 = R −d, and calculate the magnetic field at the origin. We can see that this problem is symmetric to the problem solved in the previous section, by switching between the dipole and the point at which we measure the field. Therefore, in this case we get the exact same behaviour as in Eq. 6. We can now combine both effects: shifting the dipole to distance d from one edge, and calculating the field at a distance d from the second edge, leading to a √ R improvement from each effect, resulting in Eq. 7.

NUMERICAL SIMULATIONS
We augment our analytical analysis with numerical simulations, to address an optimized aperture geometry which is not circular (but rather elliptical) and thus cannot be addressed analytically. The simulations were performed by solving the London equations as described in [39]. This method is valid for a magnetic field much smaller than the upper critical field H c2 of the superconductor, and for a superconductor much larger than its coherence length ξ.
The physical quantity obtained from the simulation is the thin film current J (x, y) = dzj(x, y, z) = (J x , J y ). If the thickness is nearly constant with thickness ∆, and the film is thin enough such that j(x, y, z) is not dependent on z, we can approximate J (x, y) = j(x, y, z)∆.
Since the divergence of the current is zero ∇ · J = 0, we can express it in terms of a scalar potential called the stream function g(x, y) J = −ẑ × ∇g = ∇ × (ẑg) = (∂g/∂y, −∂g/∂x).
The simulations find the current stream function g(x, y) and the effective magnetic field H z (x, y) resulting from an arbitrarily shaped superconducting film with an applied magnetic field H a (x, y). The simulations are carried out on a non-equidistant grid and use a matrix inversion method with matrix size (N x ×N y )×(N x ×N y ), where N x,y are the numbers of grid points in directions (x, y). This process is very resource intensive and therefore only works for small grids (usually grid sizes of < ∼ 100 × 100 points).
The simulation invert a matrix whose non-diagonal terms are multiplied by the effective penetration depth (also known as the 2D screening length or Pearl length) , where ∆ is again the superconducting film thickness and λ is the London penetration depth. Therefore, in order for the matrix inversion to be smooth, we require Λ 0. The London penetration depth λ can be smaller than the film thickness ∆ as long as it is of the same order of magnitude.
For our simulations, we choose ∆ = 80 nm and λ = 50 nm, which correspond to the London penetration depth for a clean Niobium layer of the given thickness. The superconducting film is taken to be 90 times the radius of the aperture and the entire simulation grid size is taken to be 100 times the radius. The applied magnetic field is taken to be the field created from a magnetic dipole pointing in theẑ direction: H a (x, y) = m 2πr 3ẑ . The magnetization is of a single spin 1 particle: m = 2gµ B where g is the electron g-factor and µ B is the Bohr magneton.
In the symmetric case, for which the dipole is at the center of a round aperture, there are 100 grid points per axis. The points are chosen such that there is a higher density of points next to the edges, where the gradient is higher. For the other cases, the number of points per axis is not the same. The simulations were executed with varying numbers of grid points to ensure stability and convergence (For a detailed explanation of the simulation and the choice of grid points see [39]).

Comparison between numerical simulations and analytical calculation
We first simulated the case of a round aperture [ Fig. 2(a,b)] in order to compare the numerical results with the results of the analytical analysis. As can be seen in Fig. 6(a,b), the simulations and analytical results exhibit the same trend, with small deviations. The London penetration depth λ has an effect of smearing the currents: a smaller value allows the currents to "change" on a shorter length-scale. This results in having the effective magnetic field more localized toward the edge of the aperture. This can partly explain the small discrepancy between the analytical results (which require λ = 0) and the numerical method which require a λ 0. Next, the simulations were carried out on circular apertures with varying radii, with the dipole either located at the center or at a distance of d = 100 nm from the left edge of the aperture. The magnetic field at a constant distance of d from the right edge of the aperture was extracted from the data. The resulting data can be seen in Fig. 6(c,d). The simulation results deviate from the analytics in certain cases, for which these results indicate a higher magnetic field inside the aperture. The cause of this discrepancy requires further study which is beyond the scope of this work. Nevertheless, we smooth the results by averaging the data using a moving window method, and the standard error was extracted using a moving std window.
We then fit the data using a weighted least squares regression and obtain a consistent power-law scaling to the one predicated by the analytical analysis (Eq. 6,7) and shown in Fig. 5.

Dipole at the side of an ellipse
We proceed to simulate an elliptical aperture with the x radius being a = 1000 nm and the y radius b = 100 nm. The dipole was located at a distance of d = 100 nm from the left edge. The results are shown in Fig. 7(a). It can be seen that the behavior is similar to the round aperture, but that the magnetic field is stronger. We then expanded the simulations to elliptical apertures with a constant b = 100 nm and varying a. The magnetic field at a constant distance of d from the right edge of the aperture was extracted from the data. The results [ Fig. 7(b)] were fitted to a power-law scaling, obtaining a power law of −1.4 ± 0.3.

DISCUSSION AND SUMMARY
This work addresses the challenge of achieving strong magnetic coupling between magnetic dipoles, a major obstacle for scaling spin-based qubits, and an important aspect of non-local magnetic sensing. To this end we derived an analytical solution for the flux confinement and for the magnetic field inside an aperture in a superconducting film. We have shown that for two dipoles next to the edges of a round aperture, the magnetic field, and thus the interaction, decays as 1 r 2 , which represents a significant improvement over the 1 r 3 scaling of the normal dipole-dipole interaction. We have extended these results using numerical methods to show that the scaling improves even further for an elliptical aperture, reaching 1 r 3/2 . Other structures, such as the proposed "dog-bone" [Fig. 1] is expected to exhibit similar behavior as the el-lipse due to the strong confinement in one of the axes.
A similar "dog-bone" structure, with a ferromagnet instead of a superconductor, has already been shown to improve the coupling by a ratio of L D [13], with L being the distance between the qubits and D being the radius of the "dog-bone" circular aperture. This leads to a powerlaw dependence of -2, which we have shown that can be achieved using a simpler circular aperture inside a superconducting film. By creating an ellipse or a "dog-bone" structure in a superconductor as proposed here, we expect to improve the power-law dependence even further to around -1.5. We also stress that the ferromagnetic coupler proposed in [13] introduces significant spin noise which can reduce the coherence time of the qubits, while the superconducting structure can avoid this adverse effect.
Two other works have shown direct NV-NV coupling between two NVs at close proximity. In one, two NVs at a distance of 10 nm were shown to have a coupling of ∼ 40 kHz [5]. In the other, two NVs at a distance of 22 nm were shown to have a coupling of ∼ 4.9 kHz. Using the method proposed in this work, we can achieve a coupling of ∼ 10 kHz at a distance of 300 nm. Such a coupling strength is significant compared to previous results mentioned above, and compared to NV coherence times achieved at cryogenic temperatures using dynamical decoupling schemes, reaching nearly 1 second [41]. Thus, the proposed structure could allow high-fidelity two-qubit gates (at a rate of over 1000 operations within the coherence time), while maintaining a spatial separation between qubits of ∼ 300 nm, well beyond the optical diffraction limit, allowing straightforward addressing of individual NV centers.
These results could significantly impact several fields, including SQUIDs [39], magnetic levitation [42], quantum sensing and quantum computation [3,43]. Specifically, our findings could potentially allow coherent long range dipole-dipole interactions, which could pave the way toward scalable quantum computing in solid-state spin qubit architectures.
While the experimental realization of the proposed system is not trivial, our scheme can be scaled to larger magnetic dipole and aperture sizes, to facilitate simpler demonstrations. It is possible to experimentally verify our results by scaling all of the relevant sizes to µm or even mm scales. This could be done, for example, by placing a micrometer-sized ferromagnetic particle near one edge of a larger aperture, and measuring the local magnetic field near the opposite edge of the aperture, as a function of aperture size. These experiments, which are still not trivial and are currently being pursued, are expected to be presented in future publications.
We note that while our analysis is semi-classical, it is sufficient to prove that such long-range entanglement is possible. In accordance with the correspondence principle, it is impossible for the classical field to scale differ-ently than the quantum interaction strength.
We also note that in the context of quantum processing, an open question remains regrading the coherence of the enhanced dipole-dipole interaction mediated by the superconductor. Due to the phase preserving nature of the supercurrent, we expect that barring additional noise sources, the interaction will be coherent. Moreover, the superconducting gap for an 80 nm thick Nb based superconductor is in the 300 GHz range [44], which is two orders of magnitude higher than the relevant energy scale for solid-state spin qubits (such as NV centers). This implies that excitations in the superconductor should not decohere the quantum spins. Other sources of noise arising from an imperfect superconducting layer will also affect the qubits. However, this noise will not be different from the noise affecting other superconducting circuits and should not change the scaling of the interaction. We therefore have strong indications that coherent coupling and entanglement will be realistically possible using this approach. In addition, improved fabrication techniques should be able to suppress these noise sources even further. Nevertheless, proof that coherent coupling is possible using this method will require further study and will be the subject of future work. Where, With R being the radius of the aperture.
Calculating the magnetic field for a circular aperture with a centered dipole Plugging in the expression for the current density J (r) = −m × ∇δ(r) in Eq. 2: The left term goes away due to G(r, r ) vanishing on σ, according to the boundary conditions. Now we have to calculate ∇ r G(r, r )| r =0 . From symmetry it is clear that we can't have aφ component. Our approach will be assuming that z, z are positive, and then since the vector potential is continuous we can plug in z = 0. As we also know that it would be symmetrical for negative values, we will get the solution for all z. With these assumptions, | r =0 → −1, and its z derivative will include a delta function that vanishes for z > 0. So from now on we can plug in = −1. The derivative calculation in respect to variable x yields: Now we will calculate F ± | r =0 ,D ± | r =0 , and their derivatives with respect to ρ and z : We choose φ = φ, which will give us laterρ =ρ. Plugging in F ± | r =0 , D ± | r =0 we get: And now writing the full gradient, and plugging in the derivatives of the terms: tan −1 (α(r)) + α(r) 1 + α(r) 2 ρρ + zẑ ≡ n(r) Where, Plugging in back the terms calculated in this expression and choosing orientation for m and angle φ, this reduces to: Numerical simulations Ampere's law for the field around a superconducting film can be written as: H z (r) = H a (r) + S d 2 r Q(r, r )g(r )) With H z (r) being the resulting perpendicular magnetic field, H a (r) being the applied perpendicular magnetic field, g(r ) is the current stream function and Q(r, r ) is a kernel that represents the perpendicular magnetic field created at point r by a magnetic dipole of unit strength andẑ direction at position r . The second London equation in 2D: H z (x, y) = −Λ ∇ × J (x, y) ẑ = Λ∇ 2 g(x, y) Eliminating Hz from Eq. 11 using Eq. 10 we can get (in discretized form): w j is the weight of point (r j ) which represent its 2d volume. The matrix ∇ 2 ij computes the 2D Laplacian at point r i from the stream function g(r j ) and its four nearest neighbors. By inverting this equation, we can directly extract the stream function from the applied magnetic field: With the inverse matrix K ij : Note that the dipole field is a factor of 2 higher along it's axis. When the dipole is in theẑ direction, the field vanishes on the superconductor, since no field lines can cross it. When the dipole is in thex orŷ directions, the field does not vanish since the superconductor has zero thickness and can't block the field lines. For theŷ case, the field at x > R coincide with the free dipole field and is not visible in the figure. Taking the limit of R → ∞ (or ρ → 0), we get the behavior of a free dipole, with Eq. 9 approaching the free dipole case. We see that when the dipole is in theẑ orŷ directions the magnetic field behave quiet similar, approaching infinity as we get closer to the edge, with the dipole in theẑ resulting in a bit stronger magnetic field. A dipole in thex direction results in a weaker field compared to the free one.
This matrix inversion is the most resource intensive process of the simulation and is the limiting factor of the grid size. The simulation works by first calculating the stream function g using Eq. 13 and then extracting the resulting magnetic field using Eq. 10. The stream function g(x, y) has several useful properties: