Wentao Qi ,Alexandr I Zenchuk ,Asutosh Kumar and Junde Wu,*

1 School of Mathematical Sciences,Zhejiang University,Hangzhou 310027,China

2 Federal Research Center of Problems of Chemical Physics and Medicinal Chemistry RAS,Chernogolovka,Russia

3 P.G.Department of Physics,Gaya College,Magadh University,Rampur,Gaya 823001,India

4 Harish-Chandra Research Institute,HBNI,Chhatnag Road,Jhunsi,Prayagraj 211019,India

5 Vaidic and Modern Physics Research Centre,Bhagal Bhim,Bhinmal,Jalore 343029,India

Abstract Fundamental matrix operations and solving linear systems of equations are ubiquitous in scientific investigations.Using the‘sender-receiver’model,we propose quantum algorithms for matrix operations such as matrix-vector product,matrix-matrix product,the sum of two matrices,and the calculation of determinant and inverse matrix.We encode the matrix entries into the probability amplitudes of the pure initial states of senders.After applying proper unitary transformation to the complete quantum system,the desired result can be found in certain blocks of the receiver’s density matrix.These quantum protocols can be used as subroutines in other quantum schemes.Furthermore,we present an alternative quantum algorithm for solving linear systems of equations.

Keywords: matrix operation,systems of linear equations,‘sender-receiver’ quantum computation model,quantum algorithm

1.Introduction

Feynman and Deutsch have pointed out that it would be impossible to accurately and efficiently simulate quantum mechanical systems on a classical computer.They also conceived the idea that machines ‘built of quantum mechanical elements which obey quantum mechanical laws’[1,2]might be able to process information fundamentally more efficiently than typical classical computers.Such computational power would then have applications to a multitude of problems within and outside quantum mechanics,including information theory,cryptography,mathematics and statistics.Foreseeing such a situation,arose the need of building ‘quantum computers’ and devising ‘quantum algorithms’ that could be run on quantum computers.Over the past three decades there has been a substantial growth in research in both theoretical and experimental aspects of quantum computing.Quantum algorithms exploit principles of quantum physics,and are known to offer significant advantages over classical counterparts for a number of problems.The first evidence of quantum advantage was the Deutsch algorithm [2],a simple quantum algorithm to quickly verify if a Boolean function is constant or balanced.Henceforth,quantum algorithms were paid much attention and had been significantly developed.Namely Simon's algorithm [3,4],Shor’s factoring algorithm[5,6] based on the quantum Fourier transform [7–9],and Grover’s algorithm [10] are much celebrated.

Figure 1.The Sender(S)-Receiver(R)model for various matrix operations:(a)General schematic.(b)Matrix-vector product.S1 has m rows of k-qubits and S2 has only one k-qubit row.R has(m +1) qubits where the first m qubits are part of S1 and the last one belongs to S2.(c)Matrix-matrix product.S1 and S2 have m and n rows of k-qubits,respectively.R has(m +n) qubits where the first m qubits belong to S1 and the remaining n qubits belong to S2.(d)Sum of two matrices.Both S1 and S2 have m rows.The first row of each sender has n+1 qubits,and the remaining m-1 rows have n qubits.R includes(m +n) qubits where the first m qubits belong to the last column of S1 and the rest n qubits belong to the last row of S2.The columns in S1 are enumerated from left to right,while the columns in S2 are enumerated from right to left.(e) Determinant of an n×n square matrix.There are n senders,and each sender Si is a single n-qubit row.R includes the last qubit of every Si.(f)Inverse of a non-degenerate n×n square matrix.There are n senders,and each sender Si is an(n+1)-qubit row.R includes the last two qubits of every Si.

It is possible to reduce many scientific problems for matrix formalism and for the problem of solving linear systems of equations.Fundamental matrix operations and solving linear systems of equations frequently arise on their own,as well as subroutines in more complex systems,and are ubiquitous in all realms of science and engineering including machine learning and optimization.However,with huge data sets and large dimensions of physical systems,such tasks can be practically intractable for classical computers in terms of data processing speed,storage and time consumption.A quantum algorithm from Harrow,Hassidim and Lloyd(HHL)for solving linear systems of equations [11] was devised recently,which provides an exponential speedup over the best known classical algorithms,and has led to further remarkable developments [12–14].Moreover,the HHL algorithm has been experimentally realized for some simple instances[15–18].An alternative protocol for solving systems of linear algebraic equations with particular realization on superconducting the quantum processor of the IBM Quantum Experience was proposed in [19].After the advent of the HHL algorithm,more efficient quantum algorithms and many substantial results for matrix operations have been obtained[20–24].In particular,Zhaoetal[21]have presented a novel method for performing elementary linear algebraic operations on a quantum computer for complex matrices.To simulate these matrix operations,they constructed a Quantum Matrix Algebra Toolbox,a set of unitary matrices,and used the Trotter product formula together with a phase estimation circuit.These results,added with quantum computation,can be used for solving many problems in machine learning and data processing.The scalar product of arbitrary vectors has also been investigated[21,22]using an ancilla and Hadamard operator.Stolze and Zenchuk,using a ‘sender-receiver’model,proposed an alternative scheme for the scalar product via a two-terminal quantum transmission line [25].They encoded the given vectors into the pure initial states of two different senders and obtained the result,as an element of the two-qubit receiver’s density matrix,after time evolution and a proper unitary transformation.In the case that the model does not consider time evolution,design of the unitary transformation becomes very significant for the whole quantum computation.

In this paper,based on the ‘sender-receiver’ model,we present quantum algorithms for matrix operations such as the matrix-vector product,the matrix-matrix product,the sum of two matrices,and the calculation of determinant and inverse of a matrix.We stress here that the protocols proposed in this paper do not use Trotterization which is a basic tool of algebraic calculations in [21] and is an operation-consuming process requiring multiple repetition of a set of certain operators.Instead,for each matrix operation,we prepare the unitary transformation of a quantum system which is universal for a given matrix dimension and does not depend on particular matrices to be operated.We recall that known quantum algorithms of matrix inversion [26,27] include classical arithmetics for inverting eigenvalues and therefore,they are not completely quantum algorithms.Here,we propose a truly quantum protocol for matrix inversion based on its definition in terms of algebraic complements.This is to avoid inverting eigenvalues of a matrix.Further to this,using these quantum algorithms,we propose an alternative quantum scheme for solving linear systems of equations.The proposed quantum algorithms are interesting quantum schemes for matrix operations,especially for solving linear systems of equations.In each case,the unitary operationWthat acts on the initial quantum state ρ(0) of the whole system,i.e.ρ=Wρ(0)W?is physically realizable.For instance,the multiple-quantum NMR technology can be used to achieve these operations.Regarding the running time,except for the calculation of the determinant and inverse of a matrix,the unitary transformationW,according to the Solovay–Kitaev theorem[9,28],can be approximated by the standard set of universal gates in polynomial time.

The paper is organized as follows.In section 2,we describe the sender-receiver model which serves as a general structural scheme for various matrix operations.In section 3,the detailed protocols for matrix operations including matrix-vector product,matrix-matrix product,matrix-matrix sum,calculation of matrix determinant and matrix inverse and solving systems of linear equations are proposed.A brief complexity analysis is presented in section 4,and the conclusion follows in section 5.Some additional calculations are presented in appendix A.

2.Sender-receiver model

In this model,theN-node quantum systemSconsists of two or more subsystems called senders.We denote theith sender bySi,so that.A node is actually a qubit andis the number of qubits with the senderSi,so thatThe receiver(denoted byR)doesn’t involve any additional node into the systemS;it consists of some nodes of the senders as shown schematically by the big dashed circle in figure 1(a).We denote theith sender without the receiver’s nodes byS′iand join allS′iinto the subsystemWe also employ the multiindex notation [29] as follows.LetXbe a string ofN-entries 0 and 1 indicating,respectively,ground or excited state of the corresponding spin.Let|X|be the sum of the entries ofXwhich equals to the number of excited spinsX(we call it the excitation number).In all our protocols,each sender carries a single excitation,i.e.the excitation number equals to the number of sendersn.The multi-indexMXis assotiated with the appropriate subsystemXofN(X)nodes.We conceive different ‘senderreceiver’ models for different matrix operations (see figure 1).

For the matrix-vector product,the matrix-matrix product,and the sum of two matrices,all entries of a particular matrix(vector) are allocated to a single sender.However,for the determinant and inverse of a square matrix,we allocate elements of each row of the matrix to separate senders.Thus,there arensenders for ann×nmatrix.Moreover,for calculating the inverse matrix,nauxiliary qubits are added to the model which renders assistance for reservingn2algebraic complements.

In our treatment,matrices and vectors are considered in the complex number field.The initial states of senders are all single-excitation pure normalized states|ψi〉,and the elements of matrixA(a vector is also a matrix) in a particular matrix operation are encoded as the probability amplitudes in those pure states.Hence,(see6If a matrix M does not satisfy the condition ＜1,it can be multiplied by a small enough number on a classical computer so that the Frobenius norm of the resulting matrix is less than one.Without the loss of generality,all matrices are constrained by this requirement.),whereis the Frobenius norm.The initial state of the whole system is

Note that it is the structure of the initial state that provides us all the necessary multiplications of the matrix elements placing the products into the appropriate position in the density matrix of a quantum state.All we need to get the final result is to calculate the sum of the certain elements of ρ(0) (maybe multiplied by ±1 in calculating the determinant and inverse matrix) and move the result into the desired position in the receiver’s density matrix.This can be reached after evolving ρ(0) by an appropriate unitary transformW,ρ:=Wρ(0)W?,with subsequent tracing overS′ :The unitary transformWsatisfies the commutation relation

Figure 2.The schemes for probabilistic realization of matrix operations.(a) Probabilistic method of calculating the amplitudes Ai.(b)Obtaining the output state of the receiver with Ai encoded into the (part of) n-excitation states,equation (7).

This is the general formula for the elements of the (-n)-order coherence matrix of the receiver.Our aim is to find the appropriate unitary transform for every matrix operation mentioned earlier,and present the corresponding receiver’s state in a convenient form.

2.1.Pure-state analysis

The density-matrix formulation is useful for theoretical study of the matrix-operation protocols.For practical computation,however,we turn to pure-state analysis.Let |ψi〉 be the pure initial state of theith sender,i=1,…,n.Then,after applyingW,we can write the pure state of the whole system as follows:

Now we can adopt one of the following two strategies.

1.Strategy of amplitudes.It is the probabilistic method of calculating amplitudesAi,see figure 2(a).Performing measurements on the receiver qubits with outputswe obtain the quantities |Ai|2as the probabilities of measuring

2.Strategy of state.It is the probabilistic method of obtaining the output state of the receiverRwith the amplitudesAiencoded into then-excitation receiver state,see figure 2(b).Performing measurements on qubits of the subsystemS′ with output∣0′ 〉S′results in the following state of the receiverR:

where only the first term on the right-hand side of the equation carries the useful information,while the second term has the structure

The constant γ in equation (7) provides the normalization

The state |Ψ〉Rin equation (7) can be used as an input state for other algorithms.

3.Matrix operations

In this section,different strategies are adopted for different matrix operations to design the unitary transformations.

3.1.Matrix-vector product

We encode the elements of matrixA=(aij)m×kand vectorv=(vj)k×1as the probability amplitudes into the pure normalized states

Here 0Xis the multi-index of zeros related to the subsystemX.Since the receiverRis an(m+1)-qubit subsystem,its(-2)-order coherence matrix ρ(R;-2)haselements(two excitations out ofm+1 receiver spins).But we only consider the matrix elements with indexwith 1 at theith and(m+1)th positions.Equation(12)for these elements yields

where for the sake of notational convenience the 0-parts inWare omitted here,and also in later discussions unless stated otherwise.That isLet

Pure-state analysis of matrix-vector product.–Here,we provide the pure-state analysis as discussed in section 2.Withand the initial state |ψ1〉?|ψ2〉,equation (5) takes the form

To register the result via the probabilistic method,figure 2(a),we perform measurements on the receiver qubits with the outputand obtain the quantitywhich is proportional to the square of the absolute value of theith element ofAv.

To produce the output state of the receiverRwith the vectorAv encoded into the 2-excitation state,figure 2(b),we perform the measurements on the qubits of the subsystemS′with the output∣0′ 〉S′.Then the receiver’s state equation (7)takes the form:

where γ is the normalization constant determined in equation (9).

3.2.Matrix-matrix product

We consider the matrix-matrix productABof two matricesA=(aij)m×kandB=(bij)k×s,see figure 1(c).The initial states of the two senders are

Pure-state analysis of the matrix-matrix product.–WithDΦ=ms,and the initial state|ψ1〉?|ψ2〉,equation (5) can be written as

To register the elements of the matrixABvia the probabilistic method,figure 2(a),we perform the measurements on the receiver qubits with the outputto obtain the quantitywhich is proportional to the square of the absolute value of the (il)th element ofAB.

To obtain the output state of the receiverRwith the matrixABencoded into the 2-excitation state,figure 2(b),we perform the measurements on the qubits of the subsystemS′with the output∣0〉S′.Then equation (7) yields the following state of the receiverR:

where γ is the normalization constant determined in equation (9).

3.3.Sum of matrices

Consider twom×smatricesCandD.The first row of each sender hass+1 elements,while the remaining rows havenelements.The initial states of two subsystems read

wherec00≠0,d00≠0,λ ≠0 and the extra node in the first row of each sender is marked by 0.Herek=1,2.The receiver consists of the last column ofS1and the last row ofS2as seen in figure 1(d).

There aremssuch elements which coincide with the total number of elements in either matrixCorD.Therefore,the elementscan be used to store the sumC+D.Let

is proportional to the sum of two proper elements ofCandD,whereare extra elements of the 1st row ofS1andS2,respectively,from the figure 1(d).Hence,

Pure-state analysis of sum of two matrices.–Equation (5),withDΦ=ms,Aij=λθ2(aij+bij) and the pure initial state |ψ1〉?|ψ2〉,yields

Using the probabilistic method for calculating elements of the matrixC+D,figure 2(a),we perform measurements on the receiver qubits with the outputand obtain the quantity(θ2λ)2∣aij+bij∣2,which is proportional to the square of the absolute value of the (ij) th element ofC+D.

To obtain the output state of the receiverRwith the matrixC+Dencoded into the 2-excitation state,figure 2(b),we perform measurements on the qubits of the subsystemS′with output∣0′ 〉S′.Equation (7) yields the following state of the receiverR:

where γ is the normalization constant determined in equation (9).

3.4.Determinant

To compute the determinant of a square matrixE=(eij)n×n,we encode everyith row ofEinto the pure initial state of the sendersSi,i=1,…,n,

Here,contrary to the previous investigations,there arensenders and one receiver.Consequently,a slight modification in the treatment is in order.

We only considerNR={11…1} andMR={00…0},so thatis exactly the unique element of the receiver’s (-n)-order coherence matrix ρ(R;-n)(in the rest of this paper we adopt notationXi≡XSi,for any indexX),

Pure-state analysis of determinant.–Equation (5),withDΦ=1 and the pure initial state |ψ1〉?…?|ψn〉,reads

Using the probabilistic method for calculating the determinant of the matrixE,figure 2(a),we perform measurements on the receiver qubits with the output |Φ(n)〉Rand obtain the quantitywhich is proportional to the square of the determinant ofE.

To produce the output state of the receiverRwith the det (E)encoded into then-excitation state,figure 2(b),we perform measurements on the qubits of the subsystemS′ with the output∣0′ 〉S′.Equation(7)yields the following state of the receiverR:

where γ is the normalization constant determined in equation (9).

3.5.Inverse of matrix

Like the determinant,to compute the inverse of a nondegenerate square matrixE=(eij)n×n,we encode everyith row ofEinto the pure state of senderSi.Moreover,unlike previous ‘sender-receiver’ models,we add the auxiliary module calledAuxto the senders (the last qubits of the senders).These auxiliary qubits appear inRwhich includes the two last spins of each sender.See figure 1(f).The normalized initial state of each senderSiis given by

We have already obtained the determinant.To find the inversewhereE*is the adjoint matrix,our goal is to find every algebraic complementEijofeij(the matrix element ofE).

We consider the (-n)-order coherence matrix with elementsand selectn2elements as follows.The receiverRcan be viewed as a two-column subsystem.Let the entries of the multi-indexNRassociated with the first column(thenth spins of each sender) be 1 except theith entry which is 0.Similarly,let the entries of the multi-indexNRassociated with the second column (i.e.the entries corresponding to the moduleAux) be 0 except thejth one which is 1.We denote such multi-index by,i,j=1,…,n.Thus,there aren2different elements in selected part of ρ(R;-n).These are used for storing the elements of the inverse matrixE-1.

We can also find the determinant ofEfrom ρ(R)by continuing to design the unitary transformationW.Let={1 010 …10}be the elementwhich corresponds to the states|0〉of all spins fromAuxand states|1〉of all spins from the first column ofR.Let the elements ofW(n)satisfy the constraint:

in addition to the constraint (35).Hence,

Pure-state analysis of matrix inverse.–Equation (5),withDΦ=n2+1 and the pure initial state |ψ1〉?…?|ψn〉,reads

where the subscript last means the set ofnth qubits of each sender and the subscript aux means the set of auxiliary qubits of each sender.

To register the inverse of the matrixEvia the probabilistic method,figure 2(a),we perform measurements on the receiver qubits with the outputto obtain either the quantitiesproportional to the square of elements ofE-1(if (ij)≠(00)) ordet2(E)proportional to the square of the determinant (if (ij)=(00)).

To form the output state of the receiverRwith the elements ofE-1encoded into then-excitation state,figure 2(b),we perform measurements on the qubits of the subsystemS′with the output∣0′ 〉S′.Then equation (7) yields the following state of the receiverR:

where γ is the normalization constant determined in equation (9).

3.6.Linear system of equations

We apply the quantum algorithms investigated in this paper to solve the linear system of equations,Ex=b,whereE=(eij)n×nis a non-degenerate matrix andb=(b1,b2,…,bn)Tis a unit vector.To obtain the solution x=E-1b,we first construct the unitary transformationVsatisfying the condition(2)by encoding b intoV,and then apply it to the receiver’s state ρ(R)of the inverse matrix operation: ξ(R)=Vρ(R)V?.Let

Pure-state analysis of linear system of equations.–Equation (5),withDΦ=nand pure initial state|ψ1〉?…?|ψn〉 after replacingWwith(IS′?V)W,reads:

Using the probabilistic method for registration of the elements of the inverse of the matrixE,figure 2(a),we perform measurements on the receiver qubits with the outputto obtain the quantities∣θ4σdet(E)xi∣2,which are proportional to |xi|2.

To obtain the output state of the receiverRwithxiencoded into then-excitation state,figure 2(b),we perform measurements on the qubits of the subsystemS′ with the output∣0′ 〉S′.Then,equation (7) yields the following state of the receiverR:

where γ is the normalization constant determined in equation (9).

For a summary of the proposed quantum algorithms and the illustrations of the quantum algorithms for the determinant and solving linear systems of equations,see Table B1 in Appendix B.

4.Complexity analysis and realization of W

In this section,we consider some features of our protocols.First,due to the tensor-product structure of the initial state of the senderS,all products between the elements of the matrices in the algorithms of the matrix-vector product and the matrix-matrix product,computation of the determinant and the inverse of matrix are obtained automatically and appear in the (-n)-order coherence matrix of the senderS(herenis the excitation number).Thus,to obtain the result,we have to perform the necessary number of addition and maybe multiplication by -1.This is realized by the blockW(n)of the unitary transformationW.The subsequent partial trace overS′ serves just to locate the result in the density matrix of the receiverR;this operation does not perform any algebraic manipulation.

The fact that multiplication of the matrix elements is performed automatically without special efforts reduces the number of algebraic operations imposed on the blockW(n)of the unitary transformationW.However,there is a significant problem in the realization of such a block.Apparently,the representation ofW(n)in terms of the two-level matrices [9]creates an obstacle in harnessing the above advantage.In fact,the number of ‘a(chǎn)ctive’ elements ofW(n)involved into the algebraic operation is predicted by two parameters: (i) the number ofn-excitation basis elements in the initial senders state,

(this number is determined by the number of senders and equalsdo not mistake it for the total number ofn-excitation basis elements which isand(ii) the number of entries of the receiver density matrix where the result of a matrix operation is finally placed.The latter fixes the numberK2of the involved rows ofW(n),while the former fixes the number of involved entries of each above row.Then the number of required two-level unitary operations [9] is

For illustration,let us consider the calculation of the inverse of ann×nmatrix.In this case,we needK2=n2rows ofW(n)corresponding tonexcited spins of the 2nqubit receiver,according to equation (35) of section 3.5.Each of these rows hasK1=n!nonzero entries which serve to combine all necessary products of the elements of the matrixEwhose inverse is to be calculated.The result is placed in the probability amplitudes of certainn-excitation receiver states.Thus,all in all,we only needn!n2exchange-operations instead ofn!n3operations of a classical algorithm.This is an advantage of a quantum algorithm.

Table 1.Parameters K1 and K2 for the matrix-vector product(Av),the matrix-matrix product (AB),the sum of matrices (C+D),the determinant (d et(E)),and the inverse of matrix (E-1).

However,then-excitation subspace of the state-space of the senderSincludesK1=(n+1)nbasis elements|i1〉?...?|in〉 (here |ik〉 meansikth qubit excited in (n+1)-qubit senderSk),and we needK2=n2rows ofW(n).Equation(50)then yieldsK=n2(2 (n+1)n-n2+1) 2.This is the number of two-level matrices required to encode such entries.The values of the other entries in those rows are not important for us.We have listed values ofK1andK2for various matrix operations in table 1.

Thus,the proposed unitary transformationW(n)is an example of a multiqubit transformation which cannot be effectively written in terms of elementary unitary transformations.Therefore,the development of special techniques for generating such transformations becomes relevant.For instance,Floquet engineering of Hamiltonians might be effective [31–33].

4.1.Unitary transformation via Hamiltonian evolution

In [19],the evolution of the 4-qubit chain underXX-Hamiltonian was used to generate the proper unitary transformation for solving a system of linear equations.We build on this idea and introduce the HamiltonianHthat conserves the excitation number in the system.TheXX-Hamiltonian can serve this purpose:

whereIαiis an α-projection of theith spin momentum,α=x,y,z,andDijare coupling constants.In the case of dipoledipole interaction in the strong magnetic field directed alongz-axis,we have

where φijis the angle between vector →and the magnetic field,γ is the gyromagnetic ratio and ? is the Planck’s constant.Then,at some fixed time instantt0,we can assign

Let the position of thejth node be determined by its coordinatesxj,yjandzj.Then,

Notice that there areN(N-1)/2 different coupling constantsDijparameterized by threeNcoordinatesxj,yjandzj.Of course,this number of free parameters is insufficient to parameterize all needed entries ofWin the protocol of each matrix algorithm of sections 3.1–3.6.Partially,this obstacle can be overcome by taking into consideration the set of additional nodes as qubits of the extended receiver[29].This allows us to increase the number of coupling constants (and consequently the number of free parameters) in the quantum system.See the construction ofW(2)in equation (53) for calculating the determinant of a 2 × 2 matrix in Appendix A and the linear system of 2 equations with 2 variables in Appendix B.

5.Conclusions

Based on the ‘sender-receiver’ model,we have constructed the quantum computational framework for calculating the matrix-vector product,the matrix-matrix product,the sum of two matrices,the determinant and the inverse of a matrix over a complex field.By encoding the information of the matrices into the pure initial states of the senders and performing the appropriate unitary transformationW,the final results are elements of the receiver’s density matrix which,in turn,can be considered as input for other quantum algorithms.The basic feature of our algorithms is that they don’t use Trotterization.In addition,the algorithm for calculating the inverse matrix does not use inversion of eigenvalues λi,i.e.transformation λi→1/λi.Therefore,the construction of the inverse matrix is completely quantum in nature.

Moreover,the unitary transformationWin each protocol is a universal transformation,i.e.it can be used to perform the appropriate operation with any matrices.Finally,as an application of the proposed quantum algorithms,in particular,the inverse of a matrix,we have presented a quantum scheme for solving linear systems of equations different from the HHL algorithm.The proposed algorithms also provide new insights to develop other quantum algorithms.Although the representation ofWin terms of elementary transformations is not effective,the development of methods for constructing appropriate evolution Hamiltonians might overcome this obstacle.The construction of unitary transformation via Hamiltonian evolution is discussed.

In future work,we intend to study quantum algorithms implementing the matrix elementary transformation,as it is a significant tool for studying fundamental matrix algebra problems,which can complete the calculation of the determinant and inverse matrix from a different standpoint within this paper.Further,declining the number of initial qubits and improving the efficiency of quantum algorithms are also valuable research problems that need to be addressed.

Acknowledgments

This project is supported by the National Natural Science Foundation of China (Grant No.12031004 and Grant No.12271474,61877054),the Fundamental Research Foundation for the Central Universities(Project No.K20210337)and the Zhejiang University Global Partnership Fund,188170+194452119/003.The work was partially funded by a state task of Russian Fundamental Investigations (State Registration No.FFSG-2024-0002).

Appendix A.Illustration for obtaining the determinant

For the square matrixE=the pure initial states of two senders can be written as

We need to design the unitary transformationWto find one entry inWρ(0)W?showing the determinant multiplied by a constant,that is,

LetW0101;1001=designed to ensure thatWis a unitary matrix.The commutation relation [W,Iz]=0 shows that theWis block-diagonal with respect to the number of excitations.Since the 4-qubit system contains no more than four excitations,the unitary transformationWincludes five blocks:W=diag(W(0),W(1),W(2),W(3),W(4)).We letW(0)=W(4)=1 andW(1)=W(3)=I4.The block of major concern here isW(2)written in the basis{permut|0?21?2〉}={|0011〉,|0101〉,…,|1100〉}:

After partial tracing the 1st and the 3rd qubits,we have

W(2)viaHamiltonian evolution

First of all,we remark that the construction ofW(2)in equation (A1) is not unique,and we look forW(2)as in equation (53).Only the second row ofW(2)is of interest.Thus,usingN=4 in equation (51) andt0=40 in equation (53),and requiring the following structure for the second row ofW(2),

wherea,b,andcare some constants,we obtain

witha=c=-b=-0.5i.Values ofDijin equation (A4) in view of equation(52)and equation(54)require the following coordinates for nodes:

We emphasize that 3-decimal accuracy in equation (A4)requires 6-decimal accuracy in equation (A5).

Appendix B.Illustration for solving the linear system of equations

Choose

The unitary transformationWis block-diagonal,W=diag(1,W(1),W(2),…,W(6)),and we choose the two-excitation blockW(2)in the basis {permut |0?41?2〉}={|000011〉,|000101〉,…,|110000〉} of the whole system as:

Also,the unitary transformationV=diag(1,V(1),V(2),V(3),V(4)) is block-diagonal whereV2in the basis{permut|0?21?2〉}={|0011〉,|0101〉,…,|1100〉} of the receiver is

With unitary transformationsWandV,we can find the solution of a linear system of equations.[We can chooseW(i)=IandV(j)=Ifori,j≠2.].