First-Principles Screening of Metal-Organic Frameworks for Entangled Photon Pair Generation

. The transmission of strong laser light in nonlinear optical materials can generate output photons sources that carry quantum entanglement in multiple degrees of freedom, making this process a fundamentally important tool in optical quantum technology. However, the availability of efficient optical crystals for entangled light generation is severely limited in terms of diversity, thus reducing the prospects for the implementation of next-generation protocols in quantum sensing, communication and computing. To overcome this, we developed and implemented a multi-scale first-principles modeling technique for the computational discovery of novel nonlinear optical devices based on metal-organic framework (MOF) materials that can efficiently generate entangled light via spontaneous parametric down-conversion(SPDC). Using collinear degenerate type-I SPDC as a case study, we computationally screen a database of 114,373 synthesized MOF materials to establish correlations between the structure and chemical composition of MOFs with the brightness and coherence properties of entangled photon pairs. We identify a subset of 49 non-centrosymmetric mono-ligand MOF crystals with high chemical and optical stability that produce entangled photon pairs with intrinsic G (2) correlation times τ c ∼ 10 − 30 fs and pair generation rates in the range 10 4 − 10 8 s − 1 mW − 1 mm − 1 at 1064 nm. Conditions for optimal type-I phase matching are given for each MOF and relationships between pair brightness, crystal band gap and optical birefringence are discussed. Correlations between the optical properties of crystals and their constituent molecular ligands are also given. Our work paves the way for the computational design of MOF-based devices for optical quantum technology.


Introduction
Entangled photon pair generation is a fundamental technique in quantum optics and quantum information science [1,2,3,4,5].Correlations between entangled photons can be exploited in various applications, such as quantum teleportation [2], quantum computing, quantum cryptography [6,7,8,9], and quantum imaging [10,11,12].Spontaneous parametric down-conversion (SPDC) is a widely used method for generating entangled photon pairs.SPDC occurs when a high intensity laser pump (p) interacts with a nonlinear crystal and, as a result of the interaction, pairs of low energy photons are spontaneously generated.The pair components are denoted signal(s) and idler(i) [13].The output photons must obey conservation relationship in terms of their angular frequencies and wave vectors (ω p = ω s + ω i , ⃗ k p = ⃗ k s + ⃗ k i ) [13,14].The SPDC process can be divided into types, degeneracy and direction.In type-I SPDC, the generated photons have the same polarization while in type-II, the output photons have orthogonal polarization.The system is collinear if the photons have the same propagation direction as the pump and non-collinear if they propagate in different directions.Figure 1 shows a schematic diagram of the SPDC process, highlighting energy and momentum conservation.SPDC (2ω → ω) and second harmonic generation (ω → 2ω) are conjugate processes in nonlinear optics which differ on how they are experimentally implemented [13].For a material to be SPDC efficient and have high performance over a wide range of bandwidths, it should have high nonlinear susceptibility, suitable phase matching, group velocity mismatch (GVM) for spectral purity [15], optical transparency (high band gap), thermal stability and be suitable for practical operation (low cost, availability, ease of use).The ability of nonlinear optical materials to be phase matchable depends on whether they belong to non-cubic crystal classes or not.Birefringence phase matching (BPM depends on the difference of refractive indices at specific wavelengths and is only possible for anisotropic crystal systems (non-cubic).The non-cubic crystal class is divided into groups: biaxial and uniaxial [13].Triclinic, monoclinic and orthorhombic crystal systems are biaxial, and trigonal, tetragonal and hexagonal are uniaxial.In terms of refractive indices, uniaxial crystals have n X = n Y ̸ = n Z and biaxial crystals have n X ̸ = n Y ̸ = n Z .For non-linear optical (NLO) applications it is important that the crystal is noncentrosymmetric (NCS) as they have non-zero second order non-linear susceptibility tensor (χ (2) ) which defines the efficiency of the process [13].Uniaxial and biaxial crystals can be further divided into crystal class and the list of all the crystal classes which are NCS is known.
The effective non-linearity, d eff , obtained by contracting the χ (2) with the field polarization of suitable birefringence phase matching condition, is an intrinsic property of materials and plays a crucial role in developing modern optical devices as well as entangled photon sources with high brightness (number of entangled photons) [16].The brightness of the signal scales as the square of d eff , linearly with the pump power, and linearly with propagation length (L) [17,18] when crystals have a random orientation (powder) and scales as L 2 for a single crystals.The effective non-linearity is affected by several material factors such as the electronic structure, chemical composition, optical properties, microstructure, temperature and pressure [19,9,20,21,22,23,24,25,26,27].
BBO, LiNbO 3 , KDP and other leading materials [28,29,30], have been experimentally realized for nonlinear optical process but have relative disadvantages such as absence of electron delocalization, and the inability to adjust their structures [31].Organic materials have been proposed to resolve these deficiencies but they lack good mechanical strength and thermal stability.Metal-organic frameworks (MOFs) are hybrid materials consisting of metal ions or clusters coordinated to organic ligands to form 1D, 2D or 3D structures, which have shown promise for nonlinear optical process as they can combine the benefits of inorganic and organic compounds [32,33].MOFs offer interesting possibilities in terms of synthesis of new structures with high flexibility and tunability, which allows for the optimization of their optical properties.The organic linkers in MOFs can be tuned to achieve a large nonlinear optical response and the metal centers can provide high degree of anharmonicity, which can significantly enhance the nonlinear optical response [19,31].MOFs offer high surface area, and are stable under various conditions and can be synthesized with a high degree of control over their size, shape and functionality.Several MOF structures have been experimentally synthesized and some of them have already shown potential to be used in nonlinear optics [34,35,36,37,38,39,40,41]. A recent publication reports the experimental demonstration of three-wave and four-wave mixing processes using millimeter-scale MOF single crystal with well-defined phase matching properties [42].The competitive optical nonlinearities of several polycrystaline MOF samples have stimulated the need to develop a robust computational methodology to screen the MOF databases to understand the structure-property relationship to identify and design potential MOF crystals that are promising.
In previous work, we developed a multi-scale methodology to compute entangled photon pair properties for uniaxial MOF crystals [16].We have also shown that the types of ligand and their arrangement can significantly affect the entangled photon pair properties [43].In this work, we extend our methodology to take into account the calculation of d eff for biaxial crystals with appropriate phase matching condition, GVM, spectral acceptance bandwidth, photon correlation (G (2) ) function and counting rate of entangled photon pairs for the collinear type-I degenerate SPDC.We have carried out high-throughput screening of MOFs with zinc as a metal and one-type of ligand to understand how different structures affect the properties of the entangled photon source.We have selected mono-ligand MOFs with d 10 metals, specifically Zn, because of its abundance within the MOF structure database.Also, Based on our results, we provide design rules which can be used to develop new materials which are highly SPDC efficient.We used VESTA [44] to manipulate and clean the original cif file, which was obtained from the CCDC.The software suite ISOTROPY [45] was used to identify the space group of MOF crystals.All the crystal structures were optimized using POB-DZVPP basis set with shrink points 3 3 and TOLINTEG 12 12 12 12 18.We also added DFT-D version 4 for the dispersion correction.We used solid-state DFT within the Coupled Perturbed Hartree-Fock/Kohn-Sham method (CPHF/KS) [46,47] with a PBE functional to compute the dynamical dielectric tensor ϵ ij (ω) and the dynamical second-order susceptibility tensor χ

MOF selection criteria and computational details:
(2) ijk (ω) of the optimized structures.We have used CRYSTAL17 [48] package to carry out all the DFT calculations.We benchmarked these methods in previous work [16].We also used available KDP and BBO data to compare as an additional benchmark.The TOLINTEG and SHRINK points sampling were also done to check the energy convergence.Also, the same DFT method was used in previous works for optical measurement of MIRO-101 single crystal [41] and to study the affect of arrangement of ligands and temperature on entangled photon pair properties [43].
Effective nonlinearity estimation: The estimation of the effective non-linearity (d eff ) of a given crystal is carried out in the orthogonal coordinate system X Y Z in which the dielectric tensor ϵ ij is diagonal.A crystal is considered uniaxial (with optic axis Z) if ϵ XX = ϵ Y Y ̸ = ϵ ZZ , and biaxial (with the optic axis in the plane XZ) if ϵ XX ̸ = ϵ Y Y ̸ = ϵ ZZ , and one of two relations being fulfilled: are the principal values of the refractive indices.For some biaxial crystals (crystal class 1, 2 and m), the dielectric matrix is not diagonalizable [13].In those cases, the trace of dielectric matrix gives the refractive indices.
Uniaxial and biaxial crystals can be further divided into positive or negative as discussed in the Supporting Information (SI) section 1.For biaxial crystals, it is possible that the principal values of refractive indices in the crystallographic frame (a, b, c) are different than the crystallophysical frame (X, Y, Z).Crystallophysical, optical and dielectric frame represent the same concept.Among many widely used crystals such as KTP, KNbO 3 , LBO [49], BiBO [50] and α-HIO 3 , the coordinate system X, Y, Z and a, b, c can coincide or misconcide.In KTP, the coordinate systems coincide and we can designate the coincidence as X, Y, Z → a, b, c whereas for LBO, the assignment X, Y, Z → a, c, b is commonly used [49].It is important to do the necessary transformation of the crystallographic frame (a, b, c) to satisfy the biaxial condition [49,51] (n Z > n Y > n X ) and then calculate the d eff for different polarization configurations after doing the necessary transformation of the χ (2) tensor.d eff was estimated using a specific configuration by contracting the full second-order polarizability tensor at the phase matching angles for different crystal structures.Figure 3 shows the flowchart for calculating d eff of any given crystal.Detailed description for calculating birefringence, phase matching angles, d eff , G (2) , and number of entangled photons can be found in SI section 1.Data and code are available on GitHub (Codes and data) for reproducibility.Although we have explored birefringence phase matching in our study, the MOF crystals used in our study can be explored for periodic pooling during crystal growth or using other techniques which can achieve quasi-phase matching and can much more efficient.Other types of phase matching can be used for efficient SPDC where materials lack birefringence and are non-linear or domain-disordered [52].GVM between signal and pump is important for the spectral purity which is important for heralded single photon sources and entangled photon sources [15].We have estimate GVM = 1/u s − 1/u p for all our MOF crystals [53,54], where u s is the group velocity of the signal and u p is the group velocity of pump.More details on calculating GVM can be found in SI section 5.
An important feature for applications of entangled photon pairs is the correlation time between the two generated photons (τ c ).This parameter can be defined as the full width at half maximum (FWHM) of the G (2) correlation function.This correlation function not only depends on the crystal properties, but also on the properties of the experimental set-up such as the bandwidth of the detector [55,53].For type-I SPDC the G (2) function can be written as where ν correspond to a small detuning frequency from the perfect phase matching, κ is the group velocity dispersion (GVD), L is the crystal length and σ is the bandwidth of the detector.We have used single crystal length of 1 mm for our calculations.We chose to normalize the G (2) correlation function because the correlation time of the two generated photons does not depend on the amplitude of the function; instead, it relies solely on the full-width at half maximum (FWHM) of the function.In that way we can directly compare plots to determine which crystal exhibits a wider correlation function, indicating a longer correlation time.
We also assume a Gaussian model for the transverse field profile of the pump and the two generated photons (signal and idler), collecting the down-converted light into single-mode fibers such that the transverse width of the two generated photons is the same σ 1 = σ 2 .The single-mode counting rate for photon pairs generated by type-I SPDC is given by [56] where σ p and σ 1 are the transverse width of the Gaussian profiles associated with the pump and signal or idler photon respectively.n 1 and n 2 are the signal and idler refractive indices, n g1 and n g2 correspond to the group indices of the signal and idler frequency.c is the velocity of the light, eon 2 is proportional to the monochromatic pump beam peak magnitude |D 0 p |, and P = c|D 0 p | 2 πσ 2 p /n 3 ϵ 0 is the amount of power delivered by a Gaussian pump beam.The phase mismatch ∆k is given by where ω 1 correspond to the frequency of the signal and idler photon and ω p is the fixed pump frequency.
We have used the pump wavelength of 532 nm for our study and MOFs whose band gap is greater than ∼2.33 eV.This is to ensure that there is no absorbance at 532 and 1064 nm.After applying the selection criteria on the band gap, we found 49 MOFs that are suitable for this study with exception of ECIWAO whose band gap is ∼2.28 eV.

Results and discussions
Birefringence (∆n = n x -n z ) plays a crucial role in achieving phase-matching conditions for nonlinear materials.The variations in ∆n 1064nm for all the crystals used in our study are depicted in Fig. 4A.The MOF crystal AQOROP has the highest ∆n and ONOCOL01 has the smallest value.Details about crystal class, band gap and birefringence are shown in SI section 2.
Figure 4B illustrates the relationship between ∆n and the band gap for different crystals.We find that ∆n decreases with increased band gap.Furthermore, Fig. 4C displays the correlation between the HOMO-LUMO gap of the ligand monomers and the band gap the crystals that contain these ligands.A strong correlation is evident, indicating that HOMO-LUMO gap influences the band gap of the crystal, including a red shift that results from interactions within the crystal.This observation serves as a valuable design criterion for crystals and provides the ability to control the band gap.The ability to control the band gap is very important, particularly in getting transparent crystals.Moreover, ∆n holds a key role in achieving phase matching over a wider bandwidth.These features are critical for the advancement of efficient nonlinear optical materials.The careful selection of ligands emerges as a powerful tool to fine-tune these optical properties for desired specifications.The observations discussed above can be explained by considering the polarizability (α).Polarizability refers to the ease with which the electron cloud of a material can be distorted by an external electric field and is directly linked to the refractive index.A material with higher polarizability, indicating a greater ability for charge separation, will have a higher refractive index and higher ∆n [57].Charge separation and delocalization of electrons are closely related to each other.Highly delocalized electrons lead to higher polarizability which affects the refractive index and ∆n.Furthermore, electron delocalization affects the HOMO-LUMO gap for ligand molecules and the crystal band gap.Higher delocalization of electrons, gives smaller band gaps [58].Therefore, the energy gap and birefringence are inversely related to each other.
Birefringence and phase-matching are related with each other.An increase in birefringence increases the phase-matching wavelength range which can significantly affect the nonlinear efficiency.It should be noted that if birefringence is too large, serious walk-off effect can occur and nonlinear efficiency can significantly decrease even if the nonlinear coefficients are higher [59,14].If the birefringence is too small, phasematching will not occur.We have calculated the phase-matching angles of all the crystals depending on whether a crystal is uniaxial or biaxial at 1064 nm (see Method section).We benchmarked our method by reproducing phase-matching angles calculations for KTP [60] and BiBO [50,61] as shown in SI section 3.
Dispersive properties of nonlinear crystals: Nonlinear crystals dispersive properties are of great importance in the femtosecond pumping regime.This is due to their effects on the pulse width, the effective interaction length of the pump pulses, and the conversion efficiency.Temporal walk-off given by GVM and temporal broadening given by GVD impose limits on the applicability of crystal length [62].We have estimated the GVM [54,63] and Group velocity dispersion (GVD) for the crystals used in our study.Figure 5A shows the GVM of all the crystals at a signal wavelength of 1064 nm.We have also included GVM of KDP, BBO, and BiBO using the experimental data [64,65,66] and are shown in Fig. 5A as KDPEXP, BBOEXP, and BiBOEXP.Figure 5B shows the GVD of all the MOFs as a function of the band gap.We observe that the GVD decreases with an increasing band gap. Figure 5C shows the spectral acceptance bandwidth of the phenomenon at a signal wavelength (λ s = 1064 nm) using the joint spectral function for crystals belonging to crystal class 2. Details about calculating GVM, GVD, and spectral acceptance bandwidth are discussed in SI section 5. We have also included the spectral acceptance bandwidth plots of different crystal classes and GVM as a function of signal wavelength using the experimental refractive indices data of KDP, BBO, and BiBO, including the calculated BiBO GVM, in SI section 5.For NLO materials to be technologically useful, it is not enough for a material to be highly efficient at a given wavelength but also should be phase matchable.The phase matching angle together with the χ (2) tensor determine the d eff of the crystal.Figure 6A shows the overlay of possible phase matching angles, shown with black lines, with the 2D map of d eff for all the θ and ϕ for a WOVPAB crystal belonging to crystal class 222.It is important to note that although some combination of θ and ϕ can give high d eff , those angles may not be accessible due to phase-matching conditions.Figure 6B shows the phase matching angles (θ m and ϕ m ) of biaxial crystals belonging to crystal class 222.We observe that although the crystals belong to the same crystal class, they have very different phase-matching angular profiles.The effective non-linearity is an intrinsic property of a materials and plays a crucial role in developing modern optical devices as well as entangled photon sources with high brightness (number of entangled photons).The source brightness scales with d 2 eff , and d eff is used to parameterize the effective Hamiltonian for the SPDC process.In order to calculate d eff at 1064 nm for different crystals, we have divided MOFs into uniaxial and biaxial classes.We have used the phase matching condition to find phase matching angles, θ m for uniaxial, and θ m and ϕ m for biaxial at 1064 nm, and used the appropriate polarization configuration to find the d eff (see Method Section and SI Section 2 for transformation between crystallophysical and crystallographic frame to satisfy biaxial condition).The calculated χ (2) for all the MOFs used in our study are given in SI Section 4. Figure 6C  Figure 6D shows the absolute maximum value of d eff of all the crystal as a function of their band gaps.We observe that d eff decreases with increasing band gap.This is because materials with smaller band gaps enable greater degree of overlap between the wavefunctions of electrons and holes and leads to a higher probability on nonlinear optical interactions which can significantly affect the efficiency of the non-linear process.Materials with smaller band gaps lead to higher nonlinearities as compared to higher band gap because their electronic structures enable a greater degree of optical mixing and frequency conversion.We also observe that MOF crystals belonging to crystal class 422 (uniaxial crystal) have d eff = 0 pm/V (see Fig. 6D) although it has small χ (2) as shown in Fig. S2.This agrees well with the literature as crystals belonging to the 422 crystal class will have vanishing d eff if Kleinman's symmetry is valid [67,68].
Entangled photon pair properties: The G (2) (t 1 − t 2 ) intensity correlation function gives the probability of detecting two-photon pairs at different time t 1 and t 2 (see Method section for more details).Figure 7A shows the schematic of coincidence setup for measuring the two-time correlation function G (2) (τ ≡ t 1 −t 2 ). Figure 7B shows G (2) (τ ) for some of the MOFs used in this study (see SI Section 4 all the MOFs G (2) ).The FWHM of G (2) (τ ) gives the pair correlation time (τ c ). Figure 7C shows τ c as function of band gap for all the crystals.We observe that τ c decreases with increasing band gap.This is due to higher and smaller degrees of overlap between the wavefunctions of electrons and holes for smaller and larger band gaps.The degree of correlation between the generated photons is relevant in applications such as quantum computing, quantum key distribution, quantum teleportation, and quantum cryptography [6,7,9,2].The calculated correlation time for all the MOF crystals are in femtosecond range which is extremely challenging to measure them using photodetectors as they have resolution in picosecond regime.The fiber optics play a key role in making correlation time larger to nanoseconds from femtoseconds [69,55].We find that AQOROP (∼ 27.55 fs), ECIWAO (∼ 24.87 fs) and MOFTIL (∼ 25.76 fs) have the largest autocorrelation times.The number of entangled photon pairs generated in SPDC is an important parameter for entanglement-based applications [53].Figure 7D shows the number of entangled photon pairs generated by all the crystals as function of their d 2 eff .We found that AQOROP (∼ 7.47 × 10 7 s −1 mW −1 mm −1 ), ECIWAO (∼ 6.99 × 10 9 s −1 mW −1 mm −1 ) and MOFTIL (∼ 1.69×10 7 s −1 mW −1 mm −1 ) have the largest photon conversion rates.Figure 8 shows the components of AQOROP and MOFTIL MOF crystals.All the results discussed above are with signal wavelength of 1064 nm.From all the crystals we have used in our study, we could not estimate d eff for QAMFUF01 and BEKVOD crystals although they have high χ (2) as shown in Figs.S3 and S4, because we were not able to find any phase matching angle at 1064 nm.In order for these crystals to be phase-matchable, we plotted the phase-matching angle as a function of wavelength.Figure 9A shows how θ m changes with signal wavelength for QAMFUF01 (1600 -2500 nm) and BEKVOD (2200 -3200 nm).We then calculated χ (2) for QAMFUF01 and BEKVOD crystals at 1600 and 2200 nm. Figure 9B shows the d eff of these two crystals at 1600 and 2200 nm.

First-Principles Screening of Metal-Organic Frameworks for Entangled Photon Pair Generation14
Having nonlinear optical materials to be operational over wider bandwidth is an important requirement for photonics application.We have shown in Table 1 how the photon conversion rate is affected as a function of pump wavelength.We chose MOFTIL because it is the best crystal for entangled photon sources with pair conversion rate of 1.7 × 10 7 s −1 mW −1 mm −1 , and birefringence of 0.42 at 1064 nm and band gap of 3.2 eV.We observe that d eff and entangled photon pairs decrease with increasing pump wavelength (decreasing photon energy).In summary, the band gap, birefringence, signal wavelength, phase matching angles, effective non-linearity (contracted form of χ (2) tensor), number of entangled photons are correlated with each other.Figure 10 shows how all of these quantities are interrelated and how these correlation can be exploited for designing new materials which will be SPDC efficient.First-Principles Screening of Metal-Organic Frameworks for Entangled Photon Pair Generation16

Conclusion
We have successfully explored MOF database for mono-ligand MOFs with Zn as a metal center for generating entangled light via collinear degenerate Type-I SPDC.All the MOFs in our database is experimentally synthesized and come from CCDC after applying MOF filter as discussed in the method section.We have systematically developed and implemented a multiscale computational technique to the selected MOFs (uniaxial or biaxial) to identify the features that are important for high brightness and large photon correlation times.We found MOFTIL crystal to be the best candidate entangled photon source with photon conversion rates of 1.7 × 10 7 pairs s −1 mW −1 mm −1 , and birefringence of 0.42 at 1064 nm and band gap of 3.2 eV.We observe that the band gap of the crystals is related to the properties such as birefringence (an important parameter for angle phase matching), effective nonlinearity and the entangled photon pair properties.Band gap can be related to the HOMO-LUMO gap of the ligand, assuming that the ligands do not interact significantly in the crystal, which can be a design criteria for preparing new MOF crystals with desired band gaps.Advanced entangled photon sources are crucial for photonic quantum technology and MOF nonlinear materials show great promise for modern nonlinear applications.

Acknowledgments
YJC thanks the University of Notre Dame for financial support through start-up funds.SP, RAF and FH are supported by ANID through ANID Fondecyt Postdoctoral 3220857, Fondecyt Regular 1221420 and Millennium Science Initiative Program ICN170 12.We thank center for research computing (CRC) at University of Notre Dame (UND) for computing resources.

Data and Codes Availability
All the data for Sellmeier curves and second-order susceptibility tensor obtained with CRYSTAL17 are available on GitHub.Codes to calculate phase-matching angles, effective nonlinearity, correlation function and number of entangled photon pairs for Type-I SPDC are available on GitHub.Link to GitHub page is GitHub-Data-Codes.Effective non-linearity:

References
The electric field component in the ordinary axis is given by E o j (ω) = (a j )E o (ω) and the electric field component in the extraordinary axis is given by E e j (ω) = (b j )E e (ω), where, Uniaxial crystal: For type-I negative uniaxial, d eff is calculated as and for a positive uniaxial crystal, we use P (ω 3 ) is the dielectric polarization of the crystal at the pump frequency A conventional e-oo type-I condition was used [1] for negative uniaxial crystal, where the entangled signal and idler are polarized along the ordinary axis (o) and the pump is polarized along the extraordinary axis (e).For a positive uniaxial crystal, a conventional o-ee type-I condition was used Biaxial crystals: To estimate the d eff for biaxial crystals, it is important that the chosen crystallographic axes a, b, c does not differ from the crystallophysical (optical) one (X, Y, Z), for which one of two relations (n [4,2,6,7].We have satisfied the relation n Z > n Y > n X and necessary transformation of tensor element is done to calculate the d eff using Eqns. 2 and 3 for negative and positive biaxial crystal.The phase matching direction which gives the maximum d eff among all possible phase matching directions, solving equation 5, is termed the optimum phase matching direction.We use the maximum value of d eff for other analysis. It is important to note that in order to calculate d eff , we have to know whether the crystal is uniaxial (positive/negative) or biaxial (positive/negative) and what are the phase matching angles (θ m for uniaxial, θ m and ϕ m for biaxial).The difference in the phase matching angle is due to one axis enabling symmetry of revolution in uniaxial crystals compared to biaxial crystals which lack that symmetry.All the discussion about a crystal being positive/negative uniaxial or positive/negative biaxial and their phase matching angle calculations are discussed below.The angles defined are according to the positive nonlinear optics frame convention of Roberts [7].
Positive/negative crystal and Birefringence: The value of birefringence in uniaxial crystal is the difference between n o (n X ) and n e (n Z ), ∆n = n o − n e .For n o > n e (n o < n e ), the material is termed a negative (positive) uniaxial crystals [8,9].For materials which are biaxial, the situation is very complex.The three refractive indices n X , n Y and n Z along X, Y and Z axis are different.The semiaxes of the ellipsoid can always be chosen in a way that either the material is termed as positive in analogy to uniaxial crystals, and if n Z -n Y < n Y -n X , the material is termed as negative in analogy to uniaxial crystals [8,9].It is important to note that through refractive indices, the crystal type (positive/negative) is wavelength (temperature) dependent.In a biaxial crystal, the birefringence is taken to be the difference between maximum (n Z ) and minimum (n X ) refractive index.
Sellmeier curves and Type-I phase matching angle: For the phase matching geometry for different crystals, we are interested in the Sellmeier curves along different axes.We have obtained the Sellmeier curves n 2 i=X,Y,Z (λ) for all the MOFs crystals from the dielectric tensor ϵ ij (ω), computed at wavelength intervals from λ min = hc/E g to λ max = 1100 nm.E g is the band gap for different MOFs crystals.We fit the Sellmeier curves for each crystal axis to the standard relation to obtain the Sellmeier coefficients.For a uniaxial crystal, the Type-I phase-matching angle, θ m , is taken from the optic axis.The phase-matching "point" in three-dimensional form symmetrical cones [10,1].We can obtain θ m for negative uniaxial crystal using sin ), and for positive uniaxial crystal using sin ) [10].To calculate the phase matching angles θ m and ϕ m in biaxial crystals, one must solve the following equation for θ and ϕ where k X = sin θ cos ϕ, k Y = sin θ sin ϕ, and k Z = cos θ. θ is the angle from the Zaxis and ϕ is the angle from the X-axis in the XY plane.n X , n Y and n Z are the three principal refractive indices.Equation 5 can be simplified to equation 6 by letting On solving equation 6, we get For biaxial type-I phase matching, we solve where the subscripts 1 and 2 stand for fundamental and second harmonic wavelengths.On solving Eq. 8, we get θ and ϕ for biaxial type-I phase matching.More details about type-I phase-matching in biaxial crystal are given by Yao and Fahlen [3].Glauber Correlation function: The detection of arrival times of all the photon pairs generated by the MOF crystal is done in a coincidence setup.The probability of detecting two-photon pairs at times t 1 and t 2 is proportional to the correlation function G (2) (t 1 , t 2 ) given by [11,12] where the electric field operators are given by with ϵ k,σ = ℏω(k, σ)/2ϵ 0 n 2 (k, ω).V is the quantization volume, index σ is summed over a two-dimensional orthogonal polarization vector.k is an index summed over all wave vectors, e k,σ is the two-dimensional unit polarization vector, α k,σ is the mode amplitude, G(ω) is a transmission function that models the detector filter, ω the frequency, ϵ 0 is the permittivity of free space.
We have assumed that this correlation depends only on a relative time τ = t 1 − t 2 .It is possible because we only consider the pump as a gaussian beam and the Hermite-Gauss basis for the generated photons.The right-hand side of Eq. ( 9) is proportional to the frequency spectrum of the light as and the light spectrum frequency is given by where ν is a small detuning frequency, and Φ(x) is known as the joint spectral function (JSA) which is typically proportional to a "sinc" function.∆k is given by ∆k where n g is the group velocity difference for the generated photons, and κ corresponds to the group velocity dispersion of the generated photon.Since we are interested in Type I degenerate SPDC, the generated photons have the same central frequency that is half the central frequency of the pumped photon.As the generated photons have the same polarization, n g is and the ∆k expression becomes On combining Eqs. 12 and 14 in Eq. 11, we obtain Glauber Correlation function as Counting Rate: The biphoton generation rate depends on the integral [13] |Φ The problem simplifies for the case in which a collimated Gaussian pump beam directed onto an isotropic rectangular crystal (positioned at the origin within a Cartesian coordinate system, with the ẑ axis aligned along the optic axis) with dimensions L x × L y × L z .Assuming our interest lies in gathering the down-converted light through single-mode fibers, only photons generated in the zeroth-order Hermite-Gaussian modes will impact the detected event rate.Consequently, G p (⃗ r), g ⃗ µ 1 , and g ⃗ µ 2 all conform to gaussian functions, leading to the transformation of |Φ(∆k z )| 2 to where . We have assumed that both detectors' filters possess identical bandwidths, denoted as σ 1 = σ 2 .Expanding the limits of integration over x and y to arbitrarily large values requires only that the transverse width of the crystal exceeds the dimensions of the Gaussian pump beam as well as those of the signal and idler modes.Additionally, in single-mode nonlinear waveguides, light confinement within a considerably smaller beam diameter is achievable without divergence.Under these assumptions, the overlap integral simplifies significantly to [13] |Φ Presently, an expression can be formulated for the overall down-conversion rate, denoted as R, using as a Gaussian pump beam into another Gaussian signal-idler mode, this expression is formulated in relation to the frequencies of the signal and idler, denoted as ω 1 and ω 2 as Where . From energy conservation we have ω p = ω 1 + ω 2 .Now we can integrate over one of the frequencies using the delta function.We integrate over the frequency ω 2 getting where the phase matching function is given by This function is analogous to the expression of equation 13, meaning that in the case of type I SPDC we need to replace ∆k for this expression.

Estimation of HOMO-LUMO gap:
To calculate the HOMO-LUMO gap for every ligand in every MOF, first, we retrieve the chemical name of the organic molecules from the chemical name of the MOFs in the Cambridge Structural Database (CSD) [14].By using conventional inorganic compound naming we were able to filter and obtain the chemical name of the ligand for Zinc based MOFs with one organic molecule as ligand and Zn 2+ as node.Using OPSIN [15] we converted these names to SMILES and posteriorly with Rdkit [16].We generated a 3D structure for every molecule and subsequently perform geometry optimization of the coordinates with Dmol3 [17] with Minessota funtional and DNP 4.4 basis set of Dmol3 and obtain the values for highest occupied molecular orbital and the lowest unoccupied molecular orbital, these values are used to estimate the HOMO-LUMO gap.A more refined approximation would be to extract the coordinates of the ligands from the crystal but a correct capping of each ligand is required to represent the correct valence.This chemical name-tostructure approach requires precise chemical naming on the databases, as misspelled or bad labeling is possible.

Biaxial and Uniaxial MOF properties
Table 1: Crystal system, crystal class, space group, band gap, coordinate transformation, birefringence and crystal type for biaxial crystals and that N is for negative and P is for positive.

Dispersive properties of nonlinear crystals
We have used the following functional form of sellmeier equation to fit the refractive indices values.
Here, λ is in nm.We used the fitted sellmeier equations to obtain refractive indices values at different wavelength and then estimate the GVM (see Eq. 23).

Group velocity mismatch (GVM):
The GVM for the degenerate case for uniaxial and biaxial MOF crystal was calculate using [18,12] where, u(Ω) are the group velocities at the frequency Ω.For uniaxial crystals, phase matching equations are iterated to calculate the phase matching angles at a given wavlength and then used Eq.23 to estimate the GVM.For biaxial crystals [3], we estimated the gvm at all possible phase matching angles (see SI section 1) at a given wavelength and then estimated the GVM using Eq.23.We picked the minimum value from all the absolute values of GVM. Figure 13 shows the GVM of KDP, BBO and BiBO using the experimental data [19,20,21] and are shown in Fig. 13A as KDP-EXP, BBO-EXP and BIBO-EXP.We have also included the calculated GVM of BiBO as BIBO-CAL in Fig. 13.

Group velocity dispersion(GVD):
The GVD for biaxial and uniaxial using GVD = λ 3 2πc 2 d 2 n dλ 2 , with appropriate signal depending on whether a given crystal is positive uniaxial/biaxial or negative uniaxial/biaxial (see SI section 1 and method section in the main text).
Phase matching function: Taking into account that the conditions of perfect phase matching are fulfilled and the crystal is uniaxial, we can calculate the spectral acceptance bandwidth of the phenomenon at a fixed pump wavelength λ p using the joint spectral function, which depends on the phase matching mismatch ∆k: where L is the length of the crystal.The phase mismatch can be calculated directly using a relation between the refractive indices of the pumped and the generated photons: ∆k = (n o (λ s ) − n eff (λ p , θ)) 4π λ p (25) where n o and n eff are the ordinary and effective refractive indices, respectively.We assume the crystal is positively uniaxial.For a negatively uniaxial crystal, the ordinary refractive index is calculated using the output wavelength (λ s ), and the effective refractive index is calculated using the pumped wavelength (λ p ) for SPDC (Spontaneous Parametric Down-Conversion).Consequently, the spectral acceptance bandwidth can be determined by estimating the full-width at half maximum (FWHM) of the normalized joint spectral function (24).

Figure 1 .
Figure 1.(A) Schematic of SPDC process with momentum and energy conservation.One photon energy ℏω 3 splits into two twin photons at energies ℏω 1 and ℏω 2 .(B) Phase matching configuration for type-I SPDC.(C) Sketch of MOF materials.

Figure 2
shows the flowchart of the MOF filter.After applying the selection criteria, we obtained a database of monoligand MOFs with Zn as a metal.All the MOFs in our database are experimentally synthesized and they come the Cambridge Crystallographic Data Center (CCDC).

Figure 2 .
Figure 2. Flowchart of the MOF selection filter.MOFs = 114373 represents total number of MOFs from CCDC database.Zn = 2243 represents total number of nondisordered non-centrosymmetric MOFs with Zinc as metal.

Figure 3 .
Figure 3. Flowchart for calculating d eff of any given crystal.n X , n Y and n Z are the refractive indices.n o and n e are the ordinary and extraordinary refractive indices.The crystal optic axis Z c sets a polar angle θ c and azimuthal angle ϕ with respect to the propagation direction k. d eff is obtained by contracting the χ(2) tensor with appropriate polarization configuration.

Figure 4 .
Figure 4. (A) Birefringence ∆n for different MOF crystals.(B) ∆n as a function of band gap for different MOF crystal.(C) Correlation plot between between the homolumo gap of the ligands and band gap the crystal which has this ligands.The black line represents the perfect correlation between the band gap and HOMO-LUMO gap.

Figure 5 .
Figure 5. (A) GVM of different crystal at pump wavelength 532 nm.Blue circles represent the GVM of all the MOFs used in our study using calculated refractive indices and red circles represent the GVM of KDP, BBO and BiBO using experimental refractive indices.(B) GVD as a function of band gap for different MOF crystals.(C) spectral acceptance bandwidth plot for crystals belonging to crystal class 2.

Figure 6 .
Figure 6.(A) 2d map of d eff as function of θ, ϕ of WOVPAB MOF crystal belonging to 222 crystal class.Black lines shows all the possible phase matching angles (θ m and ϕ m ).(B) Phase matching angles (θ m and ϕ m ) of biaxial crystals belonging to crystal class 222.(C) d eff of crystals belonging to crystal class 222 as function of phase matching angle, θ m .ϕ m is not shown for clarity.(D) d eff as function of band gap for different crystals.Crystals whose d eff values are closer to 10 −6 pm/V belong to the 422 crystal class.
shows d eff as function of θ (ϕ m is not shown for clarity) for MOF crystals belonging to crystal class 222.The sign of d eff is unimportant and it is the profile of d 2 eff that dictates the conversion efficiency along a given phase matching direction.Phase matching angles (θ m ) of all the uniaxial crystals are shown in SI Section 2 and phase matching angles (θ m and ϕ m ) of all the biaxial crystals used in our study are shown in the SI Section 3.

Figure 7 .
Figure 7. (A) Visible photons (green) are pumped into a NCS MOF single crystal with propagation length L. Entangled photon pairs at half the pump frequency (red), produce from the crystal with a small delay τ .Photon pairs are collected with a coincidence detector (CD), and the delay times are post-processed to give a two-time correlation function G (2) (τ ).A special case of type-I SPDC where the signal and idler are coming to the same detector and the difference between the arrival time of the photons represents the correlation time.(B) G (2) (τ ) of selected MOFs.(C) Entangled photon correlation time as a function of band gap for all the crystals.(D) Pair brightness as a function of d 2 eff .

Figure 9 .
Figure 9. (A) Phase matching angle (θ m ) for QAMFUF01 and BEKVOD crystals for a range of signal wavelengths.(B) Effective nonlinearity (d eff ) for for QAMFUF01 crystal at 1600 nm and BEKVOD crystal at 2200 nm.

Figure 10 .
Figure 10.Flowchart showing how different quantities are correlated with each other.

3 .
Figure 1.(A) Experimental [5] and calculated phase-matching angles for BiBO crystal belonging to crystal class 2. (B) d eff of BiBO crystal belonging to crystal as a function of phase matching angle, θ. (C) Experimental [3] and calculated phase-matching angles for KTP crystal belonging to crystal class mm2.(B) d eff of KTP crystal belonging to KTP crystal as a function of phase matching angle, ϕ.θ is not shown for clarity

Figure 2 .Figure 3 .
Figure 2. (Left) Phase-matching angle for MOF crystal belonging to crystal class 2. (Right) d eff of crystals belonging to crystal class 2 as a function of phase matching angle, θ. ϕ is not shown for clarity

Figure 4 .Figure 5 .
Figure 4. (Left) Phase-matching angle for MOF crystal belonging to crystal class 222.(Right) d ]rmef f of crystals belonging to crystal class 222 as a function of phase matching angle, θ. ϕ is not shown for clarity

Figure 13 .
Figure 13.Calculated GVM for KDP, BBO and BiBO using experimental and calculated refractive indices.

Figure 14 .
Figure 14.Spectral acceptance bandwidth for different crystal class.(Top left) Crystals belonging to crystal class 222.(Top right) Crystals belonging to crystal class mm2.(Bottom left) Crystals belonging to crystal class m. (Bottom right) Uniaxial crystals.

Table 1 :
d eff , number of entangled photon pairs and photon pair correlation time at different signal wavelength for MOFTIL crystal.

Table 2 :
Crystal system, crystal class, space group, band gap, coordinate transformation, birefringence and crystal type for uniaxial crystals.npm stands for no phase matching.